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ABSTRACT 

The formation of a large volcano loads the underlying lithospheric plate and can lead to 
lithospheric flexure and faulting. In turn, lithospheric stresses affect the stress field beneath and 
within the volcanic edifice and can influence magma transport Modeling the interaction of these 
processes is crucial to an understanding of the history of eruption characteristics and tectonic 
deformation of large volcanoes. We develop models of time-dependent stress and deformation for 
the Tharsis volcanoes on Mars. A finite element code is used that simulates viscoelastic flow in the 
mantle and elastic plate flexural behavior. We calculate stresses and displacements due to a 
volcano-shaped load emplaced on an elastic plate. Models variously incorporate growth of the 
volcanic load with time and a detachment between volcano and lithosphere. The models illustrate 
the manner in which time-dependent stresses induced by lithospheric plate flexure beneath the 
volcanic load may affect eruption histories, and the derived stress fields can be related to tectonic 
features on and surrounding Martian volcanoes. After an initial load increment is emplaced, 
flexurally-induced stresses grow with time and the principal stress directions in the volcano rotate 
as flexure proceeds. Magma is expected to propagate perpendicular to the least compressive stress 
axis. As a result of flexure, this axis rotates from a horizontal orientation to a nearly vertical one; 
thus magma propagation paths will tend to rotate from vertical to horizontal orientations. We 
suggest that at the later stages of flexure, this effect would tend to favor eruption sites on to the 
flanks of the volcano rather than the summit. Such a scenario is consistent with the 
photogeologically determined evolution of the Tharsis Montes. As a result of flexure there are 
three regions where stresses become sufficiently large to cause failure by faulting (according to the 
Mohr-Coulomb criterion): at the surface of the plate just outward of the volcano, and at the base of 
the elastic lithosphere beneath the center of the volcano, and on the upper flanks of the first 
volcanic load increment. Normal faulting is the dominant mode of failure predicted for the first 
region, consistent with circumferential graben observed around the Tharsis Montes and with the 
scarp at the base of Olympus Mons, interpreted as a large-offset, listric normal fault. Normal 
faulting, mostly radially oriented, is predicted for the second region. Under the premise that failure 
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in the sub-volcano zone may be strongly influenced by pre-existing weak zones or regional 
stresses, this feature may have a surface expression in the rifting and the bilateral symmetry of the 
Tharsis Montes. Failure in the third region is predicted to consist of thrust faulting, oriented 
mostly circumferentially on the upper flanks with an annulus of radially oriented thrust faulting 
starting about midway downslope. Concentric terraces, interpreted by some workers as thrust 
faults, on the upper flanks of Olympus Mons may correspond to the predicted circumferential 
thrust features. For volcanoes detached from the plate, predicted failure in the edifice takes the 
form of radial normal faulting at the base. Addition of an extensional stress due to regional 
topography yields a pattern of predicted faulting which closely matches that observed on the 
Tharsis Montes. This stress state is also consistent with an interpretation of the aureole deposits of 
Olympus Mons as the result of gravity sliding along a basal detachment. 

Our models also suggest explanations for the lack of strike-slip features, predicted by simple 
flexural models, around the Tharsis volcanoes. For a given load increment, the first mode of near- 
surface failure for most of the area immediately outward of the load is circumferential normal 
faulting and graben formation. As the volcano grows and the flexural response to the increasing 
load proceeds, the predicted failure mode in a portion of this annular region surrounding the 
volcano changes to strike-slip faulting. Because normal faulting has been predicted to have taken 
place earlier, however, it is likely that release of later stresses will occur by reactivation and 
growth of these normal faults and graben rather than by the formation of strike-slip faults. 
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INTRODUCTION 

The major Tharsis volcanoes of Mars (Olympus Mons and the three Th arsis Montes: Arsia 
Mons, Pavonis Mons, and Ascraeus Mons) are among the largest known volcanic structures in the 
solar system. An understanding of the formation and evolution of these structures can provide 
important constraints on the processes that built and maintained the Tharsis Rise. In addition, the 
Tharsis volcanoes may be analogues to large hot-spot volcanoes on Earth, such as those of 
Hawaii. For instance, Kilauea and Olympus Mons have very similar ratios of volcano height to 
diameter and of basal scarp height to volcano height [Borgia et al., 1990]. Thus studies of the 
evolution of large Martian volcanoes may yield insight into terrestrial volcanic processes as well. 

In this paper we utilize finite element models to evaluate the evolution in internal stress and 
deformation within and surrounding the Tharsis volcanoes, and we discuss how the time- 
dependent stress field may be relatable to the eruption characteristics of the volcanoes and to the 
formation of associated tectonic features. 

To date the investigation of the evolution of stresses in large volcanoes has taken two paths: 
models of edifice stresses alone (usually finite element models with rigid bottom boundary 
conditions), and investigations of flexural stresses in the lithospheric plate supporting the volcano. 
As examples of the first category of study, Chevallier and Verwoerd [1988] used an axisymmetric 
planar finite element code to investigate the effect of magma chamber and external pressures on 
stress and eruption histories of hot spot volcanoes, Dieterich [1988] modeled stress in volcano rift 
zones by means of a two-dimensional triangular grid of elements, and Ryan [1988] employed a 
horizontal planar finite element model of the flank of Kilauea to determine displacements and 
stresses due to dike emplacement. In an example of the second class of study, Thurber and Gripp 
[1988] applied a flexural model to constrain the tectonics of volcano flank motions, ten Brink and 
Brocher [1987] proposed a link between flexural stresses in the lithosphere and eruption history; in 
their scenario the orientation of flexural stresses beneath a given point on the volcanic chain 
changes with time as the volcanic load is emplaced and then eroded. Given that magma seeks to 
propagate along paths perpendicular to the direction of least compressive stress, magma ascent can 
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be blocked when both principal horizontal stress deviators are compressional. Thus eruption 
history depends on location along the volcanic chain or on time since volcano formation. 

Recent studies have also been carried out on edifice stresses of large Martian volcanoes. 
Thomas et al. [1990] investigated the tectonics of the flanks of Martian volcanoes by means of an 
incompressible finite-element model. Volcano self-loading and magma chamber effects were 
included, but lithospheric flexure was not. Thomas and coworkers found that stresses on the 
flanks of a large volcano are sufficient to cause circumferentially oriented thrust faulting, and they 
suggested that such thrusting produced the concentric terraces observable high on the flanks of 
Olympus Mons. Zuber and Mouginis-Mark [1992] treated the evolution of the surface of 
Olympus Mons, using a finite element model to calculate stresses in the summit caldera region 
caused by different magma chamber locations and geometries for comparison with observed 
patterns of faulting. 

In this paper we study the stress field within a volcano and the lithosphere upon which it rests 
as a unified system. With such a model formulation we can account explicitly for the interaction 
between edifice stresses (e.g., due to volcano self-loading) and flexural stresses (a result of the 
load induced by the volcano on the lithosphere). We include both growth of the volcano and 
viscoelastic deformation in the asthenosphere, so the problem is intrinsically time-dependent. The 
calculated displacements yield the subsidence history of the volcano. The orientations of principal 
stresses and their change with time provide important constraints on possible magma emplacement 
paths and eruption histories. With the computed stress fields, a failure criterion can be used to 
predict locations and modes of faulting within and near the volcano. After a short discussion of 
important geological and geophysical characteristics of the Tharsis volcanoes, we briefly describe 
the finite element procedure and the modelling assumptions. The results of the numerical 
computations are next presented, and their potential implications for the evolution of eruption 
characteristics and the formation of tectonic features are then compared with known constraints on 
the evolution of the Tharsis volcanoes. 
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CHARACTERISTICS OF MAJOR THARSIS VOLCANOES 

The four major Tharsis volcanoes are the largest of the martian shield volcanoes [e.g., Greeley 
and Spudis, 1981]. Each construct is composed of many overlapping flows and flow units 
erupted over a period of activity as much as 2-3 Gy. in duration [Tanaka, 1986]. The approximate 

heights and basal diameters are given in Table 1. 

The largest of the major Tharsis volcanoes is Olympus Mons (Figure 1). The principal 
tectonic features associated with Olympus Mons are (1) the summit caldera complex, consisting of 
a series of circular depressions with complex patterns of faulting, (2) concentric terraces seen on 
the upper slopes of the volcano, (3) the basal scarp, a nearly vertical cliff surrounding the volcano 
at a radius of about 300 km, and (4) the aureole deposits, which occupy a vast region dominantly 
to the northwest of the main shield (downslope from the Tharsis rise [U. S. Geological Survey 

1989] ). The summit caldera is likely the result of a collapse following withdrawal of magma from 
a high-level magma chamber [Mouginis-Mark, 198 1 ; Tuber and Mouginis-Mark, 1992]. As noted 
earlier Thomas et al. [1990] have interpreted the concentric terraces to be thrust faults. The basal 
scarp has been variously interpreted as a thrust fault [Morris, 1981], a listric normal fault [Francis 
and Wadge, 1983], and a fault-propagation fold over a subsurface thrust fault [Borgia et al., 

1990] . The aureole deposits are generally held to be disrupted landslide material derived from the 
slopes of the volcano [Harris, 1977; Lopes et al., 1980; Francis and Wadge, 1983]. Tanaka 
[1985] proposed that gravity sliding and spreading of the aureole deposits is enabled by a weak 
basal detachment between aureole and substrate. A layer containing about 10% interstitial or 
interbedded ice would be sufficiently weak to provide such a detachment. 

Tectonic features on the flanks of Hawaiian volcanoes on Earth may provide useful analogues 
to the Olympus Mons scarp and aureole. Acoustic backscatter images of the Hawaiian Ridge 
reveal that extensive mass wasting deposits (slumps and debris avalanches) flank all of the major 
Hawaiian shields [Moore et al, 1989]. A study of focal mechanisms of shallow earthquakes on 
the southern flanks of Kilauea [Thurber and Gripp , 1988] supports the hypothesis that a 
detachment surface separates the volcanic edifice and the older oceanic crust. By this hypothesis, 
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the volcano can oveithrust the underlying cnistal layer, sliding on a decoUement of weak oceanic 
sediments. The sUdes and slumps, driven by intrusions at the rift hones, detach from the volcano 

along large-offset listric normal faults [Upman et al., 1985]. 

Tectonic features observed on and around the Tbarsis Montes volcanoes differ somewhat from 

those of Olympus Mons. The principal tectonic features associated with these volcanoes are (1) 
circumferential graben on the lower volcano slopes and the surrounding plains, (2) linear 
embayments or rifts approximately bisecting each volcano along an axis trending N40“E and 
serving as sources for flows embaying the northeast and southwest flanks, and (3) lobe-shaped 
deposits to die west or norihwes. of the edifice [Scon and Tanaka, 1981; Scon « of., 1981a, b, c]. 
The Tharsis Montes volcanoes exhibit an approximate bilateral symmetry about a NE-SW-trending 
axis coinciding approximately with the line connecting their centers. Graben tend to occur on die 
norihwes. and southeast flanks of these volcanoes (Figures 2-5), U. perpendicular to radial 
direcdons broadly orthogonal to die line of die rifts. From photogeological study of Viking Orbiter 
images, Crumpler and Aubele [1978] proposed the following evolutionary sequence for the 
Tharsis Montes: (1) construction of the main shield, (2) outbreak of parasitic eruption centers on 
die volcano flanks along die NE-SW-trending axis, (3) subsidence of the summit and formation of 
concentric fractures and graben, and (4) formation of a bisecting rift along the NE-SW-trending 
axis, with rift eruptions leading to flooding of the summit depression and inundation of die rifted 
flanks This sequence is most advanced on Arsia Mons; Pavonis Mons has reached stage (3), and 
Ascraeus Mons stage (2). The morphology and tectonics of d* flanks of die Tharsis Montes are 
described in greater detail by Zimbleman and Ed g e« [1991], They propose that the lobate deposits 
formed by gravity sliding and were further modified by bod, effusive and pyroclastic volcanic 

activity. 

Prominent rift zones that radiate from the summit calderas are dominant structural features 
characteristic of Hawaiian shield volcanoes [Peterson and Moore, 1987], These zones are sttes of 
emplacement of magma drained from die summit magma chamber and indicate an environment of 
horizontal extension. The rift zones and landslides or slump deposits appear to be causally related, 
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as described above. The Hawaiian rift zones may be analogues to the radial rifts of the Tharsis 
Montes. 

The thickness of the elastic lithosphere beneath the Tharsis Montes volcanoes has been 
inferred from the radial distances of their circumferential graben by Comer et al. [1985] (see 

Figures 3-5). Preferred values of elastic lithosphere thickness are around 20 km. Concentric 

graben are not found around Olympus Mons, however, which led Comer and coworkers to 
conclude that the thickness of the elastic lithosphere beneath Olympus Mons must be much greater 
(> 150 km) than beneath the Tharsis Montes. Of course, the actual lithosphere does not behave 
perfectly elastically; rather its strength is limited by frictional failure at shallow depth and by ductile 
flow at greater depth. From strength envelope considerations for crustal and mantle material 
[McNutt, 1984], values for elastic plate thickness T e were converted to estimates of mechanical 
plate thickness T m and lithospheric thermal gradient by Solomon and Head [1990]. For Olympus 
Mons, flexurally induced curvature is small and T m is approximately equal to T e . The lithospheric 
thermal gradient is then less than 5 Kkm' 1 . For the Tharsis Montes, flexural curvatures are larger, 
T m exceeds T e , and mean thermal gradients (corresponding to the best fitting values of T m ) are in 
the range 10-14 K km 1 . Upper and lower bounds on T e for the Tharsis Montes from Comer et al 
[1985] allow a range of thermal gradients from 7-27 K km' 1 . 

METHOD 

We use the finite element program TECTON [Melosh and Rafesky, 1980, 1983] to model 
stresses and displacements in a large volcano and in the crust and mantle beneath and around the 
volcano. TECTON's capability for modeling a viscoelastic rtieology in the mantle allows us to 
model the time-dependent flexural response of an elastic lithosphere to each major increment of 
load. The program first calculates the elastic (i.e., "instantaneous") response to a load. The 
stresses and displacements arising from load-induced viscoelastic flow in the mantle are then 
calculated for a specified number of time steps. The Maxwell time (t M ) of the mantle, defined as 
the ratio of viscosity q to shear modulus p, is used as a convenient reference time scale (in 
subsequent discussion the term "Maxwell time" will refer to the mantle Maxwell time). Complete 
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descriptions of stresses and displacements were obtained at 5, 10, 15, 20, 30, 40, 70, 100, 200, 
400, 700 and 1000 Maxwell times after emplacement of each load increment. 

We assume that the problem is axisymmetric, with cylindrical coordinates r, 0, and z. Out-of- 
plane shear stresses and c 0z are then zero. We solve for r and z displacements and stress 
components c 00 , a^, and a n . In axial symmetry, two principal stresses are confined to the rz 
plane, and o 00 , the stress normal to this plane, is also a principal stress. An example of the finite 
element grid used for this study is shown in Figure 6. The displacement boundary conditions are 
that nodes on the side walls (r = 0 and r = r max ) are fixed in r but free to move in z, and that nodes 
on the bottom boundary (z = -820 km in the example shown) are fixed in z but free to move in r. 
The lower comers of the box are fixed in both directions. The volcano in the example has a radius 
of 200 km and is 25 km in height, the approximate present dimensions of Ascraeus Mons. In 
Figure 6 the volcano is the triangular region in the upper left-hand comer. The volcano rests on 
top of an elastic lithospheric plate of thickness T e . All elements in the volcanic edifice and the plate 
have a high viscosity appropriate to the lithosphere. These elements behave essentially elastically 
over the time scales considered here. All elements below the elastic lithosphere have a lower 
viscosity value, appropriate for an asthenosphere. These elements experience viscous relaxation 
over the time modeled. Material property values adopted in our calculations are listed in Table 2. 
Parameters such as density and Young's modulus differ for the crust and mantle. Elements above 
depth t<. are assigned crustal values for these parameters; below this depth mantle values are 
assigned. For most of the models discussed here, we have chosen T c = t c . We do not mean to 
imply that thicknesses of the Martian lithosphere and crust coincide; this choice serves only to 
simplify the models. We have performed additional calculations with T c > tc to explore the effect 
of a strong upper mantle on the flexural solutions. We take T c to be variously 20, 40, and 60 km. 
These values fall within the range obtained by Comer et al. [1985] for the Tharsis Montes and 
Olympus Mons. The outer radial boundary of the grid is taken a s r max = 1200 km. 

As a check on the use of TECTON for plate flexure calculations, we have compared the 
stresses and displacements calculated by the code to those from an analytic solution for thick plate 
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flexure in axisymmetric geometry [Comer, 1983]. We constructed a model of the plate only, 
without the volcano elements, similar to analytic models. This was done to avoid local stiffening 
of the plate due to the added thickness of the volcano near r = 0, which would change the resulting 
flexural profile. The load of the volcano was simulated by applying appropriate axisymmetric 
forces to the nodes on the top surface of the plate. The surface elevation of the plate at the final 
time step, relative to its value at the right boundary r = is shown in Figure 7a; also plotted is 
the analytic thick-plate solution. The difference between the two solutions is plotted in Figure 7b. 
The TECTON solution closely matches the analytic solution except for errors of less than 1% near 
r =0 and about 1 .5% near the flexural arch. 

We note that for cases in which the volcano is solidly attached to the lithosphere, the deflection 
given by TECTON at the final time step is less than that of the analytic solution near r = 0. Radial 
distances of the zero crossing and flexural bulge are also greater than those for the analytic 
solution. These effects result from a stiffening of the plate by the volcanic load itself. Near r = 0, 
the elements that make up the volcano take up some of the horizontal compressive flexural stresses 
that would otherwise be taken up by the top half of the plate. This effect may be important when 
analytic thin-plate flexure models are compared to topography in order to estimate elastic plate 
thicknesses, in that for welded volcanoes, T e would be overestimated. 

An important limitation of analytic models of lithospheric loading is that the load is emplaced 
instantaneously. This assumption implies a growth time t g much less than the mantle Maxwell time 
x M . For Mars, if - 10*1 Pa-s and p. ~ 10^ » Pa, x M ~ 10"> s or 300 yr. A more realistic model 
would allow the volcano to grow incrementally with time, simulating the way in which many 
discrete lava flows are superposed to form an edifice. We have implemented a modification to the 
finite element code that allows us to specify the properties of a given element as a function of time. 
The finite element grid is set up exactly as in Figure 6a, except that volcano elements are initially 
given a density of zero and a low but finite value of Young's modulus (the latter to prevent a 
singularity in the stiffness matrix used to calculate displacements). These elements are considered 
"off;" that is, for the purposes of the simulation, such elements are considered initially unoccupied 
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by volcano material. At the beginning of an appropriate time step, an element can be turned on 
by changing the values of the' density and Young's modulus to normal crustal values. Here we 
have simulated the incremental growth of a volcano by switching on elements in the volcanic 
edifice in thr ee stages (Figure 6b), at time zero, at 1000 Maxwell times, and at 2000 Maxwell 
times. We note here that our attempt to model a growing load is necessarily a coarse one. A 
typical terrestrial volcanic shield will consist of thousands of units emplaced over about 10 5 to 10 6 
years (this timescale may be 10® to 10^ years on Mars). Such units arc much lower in volume than 
those we model. Further, our increments are separated by only 1000 Maxwell times (about 30,000 
years for the parameters adopted here). This value was selected to be long enough to complete the 
flexural response to a given increment, abut is not meant to imply that units are emplaced that 
rapidly. 

In order to explore the hypothesis that volcanoes on Earth and Mars may have similar tectonic 
structures, we have used special features in the TECTON code to simulate the effect of a 
detachment between volcano and lithosphere as described by Lipman et al. [1985] and others. The 
effects of fault slip along such a detachment on the stresses and deformation can be modeled using 
the slippery node method implemented in TECTON by Melosh and Williams [1989]. We apply 
slippery nodes to all the nodes originally on the line z — 0, i.e., along the interface between the 
volcano and the lithosphere. This modification allows the volcano to thrust outward over the 
flexed lithosphere. 

The evolution of the stress field within and beneath the volcanic edifice in a model is depicted 
by means of symbols for the principal stresses within each element (Figure 8, symbols after 
Melosh and Williams [1989]). The first and second symbols denote situations where the azimuthal 
or hoop stress is the intermediate principal stress ct 2 . The third and fourth symbols represent cases 
where the hoop stress is the greatest extensional stress c 3 . The fifth and sixth symbols represent 
cases where the hoop stress is the greatest compressive stress Oj. 
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Given stress values, we use the Mote-Coulomb failure criterion to estimate regions where 
faulting would occur. The Mohr-Coulomb failure equation relates shear stress T at failure to 

normal stress a n 

Tfaiiure=c + antan<|> 


( 1 ) 


where c is the cohesive strengd, of the rock and * is tire angle of internal friction. We adopt values 
for c (3.8xl0 7 Pa) and (49°) appropriate for basalt [Handm, 1966). Elements where the shear 
areas has exceeded the Mohr-Coulomb criterion are shaded to indicate failure. Once we have 
found where failure is expected, the orientations of principal stresses are used to determine the 
style and orientation of faulting, according to the criteria of « ,1951,. Given the princtpal 
stress directions, the type of faulting (normal, thrust, or strike-slip) and the orientation (radtal or 
circumferential) can be determined. Figure 8 gives types and orientations of faulting in elements 
with the given orientations of the stress symbols. This classification scheme holds rf the pnncpal 

stresses approximately correspond to the horizontal (o,,. One) and vertical stresses ( a ) 

here If the stresses are not aligned in this way, fault geometiy may he less simply described. 
Strike-slip faults do not quite fit the usual definitions of radial and circumferential, because such 
faults are expected generally to strike obliquely to the principal stress directions. Formost 
situations, strike-slip faults make smaller angles with die o, direction than with the o 3 dtrection. 
Thus, when o, = O e9 we will consider tire resulting faults approximately circumferentially 
oriented, and when o, = <s„ we consider them approximately radial. Care must be taken when 
applying these results to regions a, depth, or where principal stiesses are oblique to the surface. 

For a general state of stress in a plane layer constrained laterally and subjected to self-load, ng 
(before die volcanic load is applied), tire orientation of tire stress symbol will be as in the firs, hne 
of Figure 8: maximum compression vettical and maximum extension horizontal. The magnitude 
of die deviatoric stress will increase with depth. This result can be derived from die equations for 
uniaxial strain (compaction) in die z-direction [e.g., Turcone end Schuben , 1982 , p. 108]: 
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<V = a 00=(^) a zz (2) 

For Poisson's ratio v = 0.25, the factor v/(l - v) is 1/3. Thus, the vertical stress will be three times 
as great as the horizontal stress, and with vertical stress equal to the overburden pressure (o^ — 
p c gz, where p c is the density of crustal material and g is the gravitational acceleration), the stress 
difference will increase with depth. The state of stress in the terrestrial crust is still a matter of 
debate. McGarr [1989] argues that a stress state in which all principal stresses are equal to the 
overburden load is a better reference stress state than the "Poisson" stress state described by 
equation (2). However, it seems likely that the state of stress lies somewhere between the two 
extremes. In general, the stress state described by equation (2) is the initial state in the lithosphere 
for most of our models. We consider additional models, however, to determine the effects of an 
isotropic lithospheric prestress (with a n = age = a a = p c gz)- For all models, the prestress in the 
asthenospheric mantle is taken to be isotropic. 

When applying the results of these calculations to actual geologic and tectonic features seen on 
Martian volcanoes, one must keep in mind the model limitations. The simplest calculations that we 
perform here start with the instantaneous placement, at time t = 0, of a significant portion of a 
volcanic edifice, equivalent to the condition t g « x M . This assumption is valid only for 
individual flows constituting at most a minor mass fraction of the volcano. While we have 
modified the finite element code to accommodate changes in element properties with time, as noted 
above the individual load increments are large fractions of the volcano mass. Other limitations are 
also noteworthy. It is difficult to determine the geometry of a given volcano load, even at the end 
of the last major shield-building eruptions, because of subsequent flexure and deformation that 
may have significantly modified the volcano's characteristics. These calculations are performed 
under the assumption of axial symmetry; possible effects of nonaxisymmetric loading (such as 
regional stress) and complex three-dimensional geometry are not incorporated. The effects of 
magma pressure, transport, and evacuation are not addressed in these models, these processes 
likely have important influences for caldera and flank tectonic evolution [Thomas et al., 1990; 
Zuber and Mouginis-Mark, 1992]. Further, lateral variations in material properties (due to 
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horizontal temperature gradient, for example) are not included. Our cbo.ee that T e - W, 
made for convenience, is not necessarily valid for the Tharsis region. Finally, it mus 
remembered tot a viscoelastic model cannot account for relief of stress due to faulung (excep 
f au ]ts modeled . prion such as in to detached volcano models). Once faulting occto, to stiess 
fields calculated ate no longer strictly valid, since faulting would relieve stresses locally antooul 
thus alter to predictions for failure at subseguen, times. As discussed below, some of 

limitations can be relaxed in future modeling efforts. 

NUMERICAL RESULTS: STRESS, DEFORMATION, AND FAULTIN 

Models with Growing Loads 

An illustrative mode, for to evolution of volcano-related stresses is shown in Figure 9, which 

depicts to case of an incrementally grown volcano on a lithosphere with T e = simple self- 
(Figure 9a) deviatoric stresses display to orientations and magnitudes expect 
loading of horizontal layers everywhere but along to top layer of volcano elements, 

teaching to upper part of to underlying plate (Figure f). 

locations and types of faulting as functions of time. For to initial volcano of 120 1cm radius, 
faulting is initiated in two areas by about 20 Maxwell times (Figure 9b). A subsurfaceregion of 
normal faulting at a depth of 35-40 km (at to base of the elastic ^ 

center of to volcanic load. Faulting also starts at to plate surface just ou “ 

th, region aU elements but one are predicted to experience circumferentially onented norm* 
faulthig (to innermost element is predicted to fail by strike-slip). By 100 Maxwell times (Figure 
9c^ both regions have expanded radially and vertically, and a new zone of faulting 
the upper volcano flanks. The pmdicted fault style and orientation on to flanks is predomman 

region. By 1000 Maxwell times (Figure 9d), to area of circumferential norma, faulung exten s up 
to 300 km radius (not shown in figure) and up to 10 km depth. The area of strike-slip faulting has 
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grown in radial extent at the expense of the region of circumferential normal faulting, but each of 
the elements in this rone has already been predicted to fail by a circumferential normal mechanism 
at an earlier time step (see Figures 9b, 9c). The upper flanks of tire volcano are predicted to fad by 
circumferential thrust faulting, and the lowermost element on the flank by radial throst faulting. 

The failure region at the bottom of the plate now extends slightly beyond the radius of the volcano 
itself. Near r = 0, this region includes the entire bottom half of tire plate (20 to 40 km below the 

original surface). 

In Figures 9e and 9f we see the effect of the next increment of volcano growth. Notice that 
most elements on tire volcano flat* predicted to fail in Figure 9d are no longer predicted to be in 
failure in Figure 9e, but some of them are predicted to fail again by the time depicted in F.gure 9f. 

No failure is predicted, however, in the layer of elements constituting tire new surface of the 
volcano. In Figure 9g, the effects of the final load increment can be seen. Only two elements ,n 

the entire edifice have deviatoric stresses sufficiently great to cause failure. 

The magnitude of the maximum deviatoric stress (the difference between the principal stresses 
of greatest and least magnitude) in the vicinity of tire volcano is shown versus time in Figure 10. A 
comparison with the stress orientation symbols of Figure 9 indicates the principal stresses used to 
calculate tile maximum deviatoric stress a. a given point. Figure 10a shows tire state of deviatonc 
stress of the elastic solution. The stress difference increases with depth, as discussed above. A.a 
given depth, tire deviatoric stress is higher beneath the initial volcano load increment than outward 
of the volcano. In Figure 10b(au= 1000 t M ) concentrations of high deviatoric stresses have 
developed in the three characteristic regions of failure from Figure 9. Deviatoric stresses exceed 
150 MPa near the volcano summit and in a wide area near the bottom of tire plate beneath the 
volcano. Failure is predicted in both areas. Deviatoric stresses are not quite so high in the near- 
surface region immediately surrounding the volcano, but faulting is predicted there because of the 
low confining pressure. There are two imporiant regions of relatively low deviatoric stress; 
directiy beneath the volcano summit near * = -4 to -14 km, and far from the volcano (at both die 
top and bottom of the plate). In Figure 10 c (at r = 2000 c„) the region of high stress in the edifice 
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in magnitude, although the vertical gradient of stress is greater for the thinner plate. Peak 
deviatoric stresses in the volcanic edifice are significantly higher for the thinner plate. 

Results for a model with a thicker elastic lithosphere and crust (T c = tc = 60 km) also reveal 
interesting changes from the first model. After the first loading increment (t = 1000 t m , Fig. 13a), 
faulting has initiated at the lithosphere base and outside the volcano, but not in the edifice itself. In 
the outer faulting region we again observe elements where predicted failure evolves from initial 
normal to later strike-slip. The area of predicted strike-slip faulting is much larger m extent than 
for the cases with thinner lithospheres. After the final load increment (Figure 13b) the region of 
predicted near-surface normal faulting expands outward to 590 km from the center (Table 3), as 
well as down to a maximum depth of 20 km. The region of failure beneath the volcano is found at 
a deeper level (30-60 km depth), because of the greater plate thickness. Faulting of any type is 

absent from the edifice throughout the period modeled. 

Maximum deviatoric stresses for this case are shown in Figure 14. Compared to the case with 
T e = 40 km (Figure 10), there are several differences. At the final time steps shown (Figures lOd 
and 14b), the region of high stress near the base of the lithosphere (as defined by, say, the 300 
MPa contour) is greater in vertical and lateral extent for the case with T e = 60 km, even though the 
maximum stress deviator is about the same. The deviatoric stress minimum, beneath the volcano, 
is also larger in extent and centered more deeply (about -14 km) than its counterpart for T e = 40 
km. Deviatoric stresses within the volcanic edifice are significantly lower for this case compared 

with the other two. 

To investigate the effects of different crustal and lithospheric thicknesses, we performed a 
calculation with a 40-km-thick crust and a 60-km-thick elastic lithosphere. The upper mantle 
lithosphere has a higher value of Young's modulus and should cause the plate to behave more 
stiffly than if the plate consisted entirely of crustal material. Figure 15 shows plots of stress 
orientation and failure for this case; compare with Figure 9 for a case with crust and lithosphere 
both 40 km thick and with Figure 13 for a case with crust and lithosphere both 60 km thick. The 
evolution of the stress field and zones of predicted failure is qualitatively similar to those discussed 
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has intensified, but remains centered at about the same position as before, at the top of the first load 
increment. The second load increment, emplaced atop the first, shows a much lower level of 
deviatoric stress than evident below. The deviatoric stress low beneath the volcano has shrunk in 
size and lowered its mean level to about z — -10 km. By the time subsidence is nearly complete 
(Figure lOd), most of the region near the volcano experiences deviatoric stress in excess of 200 
MPa. The regions of high deviatoric stress have expanded to maximum size, and the deviatoric 
stress minimum below the volcano center has again shrunk slightly. 

Changing the plate thickness can affect the evolution of regions of failure. For T c = tc = 20 
km , the entire flank is predicted to fail by thrust faulting after 1000 Tm (Figure 1 la). 
Circumferential thrust faults on the upper flanks will be surrounded by a small annulus of radially 
oriented thrusts on the lower flanks. The region of failure outward of the volcano is smaller in 
radial extent than for the case with T e = 40 km (Fig. 10b). While strike-slip faulting in the inner 
portion of this region, initial predicted failure for most of these elements occurs by circumferential 
nor mal faul tin g. In any case, this inner annulus will be covered by volcanic material later in the 
growth of the volcano. At the final stage (Figure 1 lb, t = 3000 x M ), the subsurface region of 
failure reaches through almost the entire thickness of the plate to the base of the volcano. At the 
radially distal edge of this region, stresses are rotated such that strike-slip faulting is expected, and 
the region merges with the failure region outside the volcano to form a continuous zone predicted 
to be faulting. More elements (5 vs. 2) in the edifice retain stresses large enough to cause failure 
than for the case with T e = 40 km. 

The evolution of maximum deviatoric stress for this case (Figure 12) follows a course similar 
to that for T e = 40 km. The regions of high deviatoric stress on the volcano flanks and at the 
bottom of the plate closely surround a stress minimum whose center lies at about z = -3 km. After 
the first load increment is accommodated (Figure 12a), peak stresses in the edifice are nearly twice 
as large as for the case with Tj = 40 km. A comparison of the late-stage stresses for these two 
cases (Figures lOd and 12b) shows that the deviatoric stresses near the base of the plate are similar 
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above. One effect of the stiff upper-mantle lithospheric layer is to concentrate tnto a narrower 
depth interval rite zone of failure beneath rite volcano, seemingly shielding the less ngi . 
above from failure. In Figure 1 3, this zone extends from 30 to 60 km depth, whereas in g 
only elements in rite mantle portion of rite lithosphere (40 to 60 km depth) show stresses that 
satisfy rite failure criterion. We also observe that the size of the zone of near-surface failure 
outward of rite volcano, and the amount of strike-slip faulting within it, are both greater 
plate with a strong upper mantle lithospheric layer (Table 3). Contours of — » dev 

♦i vtVir«cnhf*re Also the deviatoric stress minimum 
concentration of stress in the strong mantle lithosphere. Also, 

beneath the load occurs at a slightly shallower depth. 

Models with Instantaneous Loads 

I, is inshuctive to consider a model with an instantaneously emplaced load, both as a 
comparison with analytic models and as a reference to study die effecte that incremental V ° 
growth has on the evolution of stresses in volcanoes. Figure 17 shows stress onentauonsand 

faite eformiscaseattheconc, usionof flexure, = 1000 e M ). These are ^y 

those for dte case with incremental loading (Figure 9g) except for me volcano flanits, 

predicted to have failed along me entire surface and to 15 km depth near the summit. 
Circumferential thrust faults are predicted on the upper flanits, and an annulus of radt y 

increment of a growing volcano (Figures 9d, 1 la). The max, mum deviatonc stresses for dns case 

(Figure 18) resemble those for me incrementally loaded case (Figure lOd) except or e 

oon^tration of stresses at the top of the edifice. In this respect rite model behaves like anandyhc 

pHtTflexure model, in which deviatoric stresses are a maximum at the top and bottom 

wim a minimum at mid-plate. The volcano serves effectively to stiffen me plate near r ==0, * 

elevating the mid-plane from z = -20 km (for a 40 -km-thick plate) to z = -8 (for a plate effectively 

65 km thick at r = 0). 
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Models with Detached Loads 

The evolution of stress orientations for an incrementally grown volcano detached from the 
lithosphere is shown in Figure 19. Immediately after loading, stresses in almost the entire edifice 
are oriented with maximum compression vertical and maximum extension out of the plane (Figure 
19a). This geometry may be compared with that in Figure 9a, where both the maximum and 
minimum principal stresses are in the r-z plane in the upper flanks of the volcano, and the 
maximum compressional direction has rotated toward the horizontal. In Figure 19a the onset of 
faulting is predicted in two areas; at the plate surface just outside the volcano (as observed m most 
of the other models) and near the base of the volcano near r = 0. Radial normal faulting is 
predicted for this second region; the entire lower part of the volcano (near the detachment) exhibits 
horizontal extensional stress and vertical compressive stress, in contrast to the environment of 
generally horizontal compression in the models with welded volcanoes. This environment persists 
as the calculation proceeds (Figure 19 b-c). In Figure 19b, after 1000 Maxwell times, the area of 
radial normal faulting reaches outward through the edifice, above the detachment, to the surface at 
the volcano edge. Also note that stresses near the summit of the volcano have rotated to a state of 
horizontal compression. By the time the flexural response is complete (Figure 19c, 3000 Maxwell 
times), an annulus (with inner radius 170 km and outer radius 200 km) of radial normal faulting 
has formed on the lower flanks. Elsewhere, the stress state depicted in Figure 19c is similar to the 
final result of the welded case (Figure 9g), except that the area of predicted faulting outside the 
volcano is deeper, and the faulted area beneath the volcano does not extend as far vertically. These 
effects may be related to local stiffening of the plate by the welded volcano, which would raise the 

midplane of the plate. 

A plot of maximum deviatoric stress for this case (Figure 20) shows how the distribution of 
stresses is changed by the detachment. After first load increment, deviatoric stresses in the volcano 
tend to increase from top to bottom, the reverse of the pattern for the welded volcano (Figure 10b). 
A region of high deviatoric stress forms beneath the detachment (with approximately the same 
magnitude as the high stress region at the top of the edifice in Figure 10b). Figure 20b shows 
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stresses after completion of loading; stresses in the edifice are greatest at the bottom, directly above 
the detachment. The upper parts of the edifice have low levels of stress, except for a small area of 
slightly elevated stresses near the summit. The stress maximum beneath the detachment has grown 

in magnitude and area. 

In order to address the role of the initial stress state of the lithosphere in determining the time- 
dependent evolution of stresses and displacements, we consider a lithosphere with an isotropic 
prestress equal to the overburden pressure. For a model otherwise similar to that of Figure 9, the 
first faulting in the region outside the volcano (Figure 21a) is predicted to occur by strike-slip, in 
contrast to the prediction of circumferential normal faulting for the model of Figure 9. Notice that 
the onset of faulting in this region occurs at a later time, and the onset of faulting in the region 
below the volcano is delayed until about 400 t m (compared with 20 l M for the model in Figure 9). 
These delays are to be expected, since deviatoric stresses must build up from zero, as opposed to 
cases which start with a sizeable deviatoric stress. The extent of faulting at the final time step 
(Figure 21c) is similar (to 470 km radius, Table 3), and so is the predicted geometry of faulting in 
the region beneath the volcano (radial normal). In the early stages of flexure, however, strike-slip 
faulting is predicted for a greater proportion of the near-surface region than for the model of Figure 
9 (Table 3). After the flexure following the final load increment the proportional area of predicted 
strike-slip to circumferential normal faulting at the surface is still about 10% greater than that of the 
model of Figure 9 (Table 3). Thus we conclude that the initial stress state of the lithosphere prior 
to volcanism can affect the distribution of near-surface faulting surrounding the volcano. 

Maximum deviatoric stresses (Figure 23) far from r = 0 are small, as might be expected in a plate 
which starts out with isotropic prestress. However, the stresses near the base of the plate beneath 
the volcano and in the volcanic edifice are qualitatively similar to those for cases with Poisson 

prestress. 

DISCUSSION: IMPLICATIONS FOR VOLCANO EVOLUTION 
The time-dependent displacement and stress fields calculated in our models can be used to 
predict the type and orientation of tectonic features that would result from loading of the lithosphere 
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by a volcanic edifice. We may compare these predictions to observed tectonic features on and 
around the Tharsis volcanoes. The evolution of tectonic deformation is also constrained by the 
model results, an advantage of the finite element method over analytical flexure and other 
lithospheric loading models which consider an instantaneously emplaced load. We may then 
evaluate the validity of each model and the applicability of terrestrial analogues for some tectonic 

and volcanic features. 

Flank Tectonics: Circumferential Graben and Predicted Strike-Slip Fatdtmg 
As noted earlier, circumferential graben are observed on the lower flanks and the surrounding 
plains for all three of the Tharsis Montes (Figures 2-5). Formation of these circumferential graben 
requires an environment of radial horizontal extension at die top of the crest surrounding the 
volcano and in the near-surface region of the volcano flank. Almost every model predicts 
circumferential normal faulting in an annulus around the volcano. Models in which the volcano is 
welded to the underlying plate give near-surface horizontal compression in the edifice after 
lithospheric flexure, while models which include a detachment between volcano and plate show 
hor izontal extensi on in the lower volcano flanks. At the base of the detached volcano, the hoop 
stress dee is extensional, and the radial normal stress is congressional but small compared with 

die vertical normal stress er„. In order for circumferendal extensional features to form, <v should 
be the least compressive stress. However, since the magnitudes of a„ and <r 9e are small in this 
region, a small perturbation due to local or regional nonaxisymmetric stresses may determine the 
orientation of failure features [Nakamura, 1977]. The Tharsis Montes are located slighdy below 
die crest of die Tharsis Rise, on abroad slope to die norihwest [U. S. Geological Survey 1989], 
This long-wavelength slope of Tharsis, in the presence of a detachment, would be expected to add 
a northwest-southeast extensional stress; superposition with the axisymmetric model stress field 
would result in radial extensional features in the northeastern and southwestern quadrants (and, at 
sufficiendy high strains, a NE-SW trending rift), and circumferendal graben on the northwestern 

and southeastern flanks. 
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On Arsia Mons, graben are prominent on the half of the volcano centered on the northwestern 
azimuth (all the way up the flank to the caldera rim. Figures 2, 5), but are generally absent on the 
upper southeastem-centered half. Several graben appear at the southeastern edge of the edifice, 
however. These graben are likely to be related to extensional stresses at the surface of the 
lithosphere surrounding the volcano, as described above, and thus not related to stresses in the 
edifice per se. The graben on the northwestern half terminate near the northeast-southwest- 
trending rift that bisects the volcano. We suggest that a portion of the northwestern half of the 
volcano was able to slide down the slope of the Tharsis rise, creating a stress state that led to the 
formation of the numerous graben. In the southeastern half, however, buttressed by the rising 
slope in that direction and thus unable to move, the stress state remained stable (horizontal 
compression at the surface, as in Figure 9f). This scenario would account for the dramatic 
difference in faulting character on the two "halves" of the volcano. Of course, stresses resulting 
from caldera formation and magma chamber formation may have had a strong effect on the stress 

field near the summit as well [Zuber and Mouginis-Mark 1992]. 

We can apply the above arguments to the remaining Tharsis Montes. Graben on the flanks of 
Pavonis Mons are in general narrower than those on the northwestern flank of Arsia Mons, but 
more widely distributed in azimuth relative to Arsia Mons [Zimbleman and Edgett , 1991]. The 
widest graben are on the west and northwest flanks, analogous to the situation on the flanks of 
Arsia Mons. The graben on the north, northeast, and east flanks of Pavonis Mons, however, 
indicate that the stress state in the edifice is too complex to be described solely by the above 
scenario for Arsia Mons. Ascraeus Mons shows the least amount of normal faulting on its flanks; 
the few concentrations of graben he on the lower slopes, far from the summit caldera. Thus 
Ascraeus Mons shows the least evidence of tectonic modification due to the effects of a 
detachment Perhaps the detachment is smaller in extent there because of inherent differences in 
crustal properties, or because it has had less time to develop. The latter scenario is consistent with 
the conclusion of Crumpler and Aubele [1978] that Ascraeus Mons is at an earlier stage in its 
development than the other two volcanoes. 
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Stresses from both analytic flexure models and our models with instantaneous volcano loads 
predict strike-slip faulting at the surface of the plate, immediately outside the load Evidence for 
strike-slip faulting around any of the Tharsis volcanoes, however, is lacking. The absence of such 
features has been discussed by Golombek [1985] and Schultz and Zuber [1992], in the general 
context of large axisymmetric loads on planetary lithospheres. Examination of the plots of stress 
orientation and failure from our models (Figures 11, 13, 15) however, reveals that the predicted 
failure modes for elements in the region surrounding the volcano can change with time. In this 
region, every element that eventually fails (except the innermost) is predicted to do so first by 
concentric normal faulting. This can be seen more clearly in a Mohr diagram of the stresses in a 
typical element in this area (Figure 22). The Mohr circles grow in size as flexure proceeds, until . 
the failure envelope is exceeded at about 70 tm after first loading. The largest circle shows stresses 
at the conclusion of flexure. This stress state corresponds to that given by analytic plate models (in 
which flexure occurs instantaneously). Applying a failure criterion simply to the final stress state 
can be misleading. Consideration of the stress history in our models indicates that at the time of 
first failure, the shear failure criterion is satisfied in the surface elements immediately outside the 
volcano, and the principal stress directions predict circumferential normal faulting. These first- 
formed faults would tend to relieve extensional stresses in their vicinity and would serve as 
preferred planes of weakness for the relief of stresses accumulated later. Strike-slip features would 
thus tend to be prevented from forming. We also note that for the lower range of values of 
lithosphere thickness considered here (20-40 km), the zone of strike-slip faulting predicted from 
the late-stage stress models lies at the edge of the volcanic construct, which has the lowest surface 
elevation (Figure 7). This topographic low could accumulate lava flows and erosional deposits 
which would tend to cover deformational features (circumferential graben as well as strike-slip 
features). 

Horizontal Compression in the Edifice 

Circumferential terraces high on the flanks of Olympus have been interpreted as concentric 
thrust faults by Thomas et al. [1990], who applied elastic finite element models to investigate 
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edifice stresses and deformation. Our simplest model (a welded, instantaneously emplaced 
volcano) predicts a similar stress state, with circumferential thrust faulting on the upper flank and 
radial thrusts on the lower flank (Figure 17). Plate flexure models also predict large horizontal 
c omp ressive stresses in the area of the load, which increase with height above the mid-plane of the 
plate, reaching a m aximum at the surface, where thrust faulting is predicted. Models with growing 
loads, however, display a different distribution of stress. The maximum horizontal stress occurs 
in the lowest (first) load increment and decreases upward with each later load increment (Figures 9- 
14). Each successive load increment "feels" tectonic stresses from itself and all later load 
increments. The last units emplaced are sensitive only to the effects of their own loading and 
therefore exhibit a low degree of horizontal compression. Thus, thrust faults (predicted by simpler 
models) associated with such stresses would be expected to occur either early in the evolution of 
the volcano or not at all and, if formed, would be covered by later units which remain unfaulted 
because the failure criterion is not satisfied (Figures 9, 11). We note, however, that stresses in 
these elements are close to failure. Local stress perturbations and/or variations in internal friction 4> 
could result in some failure visible at the surface. On the other hand, if loading increments are 
smaller than those adopted here (which is likely), then the degree of horizontal compression from 
such elements would be smaller than those shown here, and faulting would not be expected. 

The presence of graben high on the flanks of the Tharsis Montes is not consistent with the 
state of horizontal compression in the edifice found in the welded volcano models. Other factors, 
such as a basal detachment, lithospheric rifting, or caldera stresses must be invoked to explain the 

presence of these faults on the upper flanks. 

Magma Transport 

The volcanic edifice stress fields calculated here have important implications for magma 
transport and eruption. Magma propagates through the lithosphere along fractures that form 
perpendicular to the direction of least compressive stress. Immediately following emplacement, the 
stresses in a given load increment are characterized by CT 3 horizontal (e.g., Figure 9a). Thus show 
that magma can propagate vertically to the summit region. As a result of flexure, however, the 
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principal stress directions in that load increment rotate such that CTi is horizontal and 03 is vertical. 
Magma propagation should tend to be horizontal in such 9 layer, leading outward from the summit 
to the flanks. Stresses in earlier load increments remain in horizontal compression as volcano 
growth proceeds. This behavior results in a maximum in horizontal stress near the core of the 
volcano, with lessening horizontal compression at higher (later) load increments. This stress 
pattern contrasts with that from analytic flexure calculations, where stresses increase linearly with 
distance from the plate midplane. Thus an implication of the models treated here is that within an 
interval of shield building there should be an evolution in favored eruption location from summit to 
lower flank. Such an evolution is consistent with the first two stages in the sequence of major 
events for the Tharsis Montes determined by Crumpler and Aubele [1978]. 

The stress evolution of a growing volcano may affect the location of magma collection zones 
or chambers. The maximum horizontal compression is achieved within the earliest load 
increments, near the bottom of the load. Such a region of high compressive stress at lower levels 
in the edifice could inhibit ascending magma from reaching the summit, and could thus lead to the 
pooling of magma. Since magma will propagate perpendicular to the direction of least compressive 
stress, radial propagation of magma in sheet dikes or sills might also occur. On the other hand, 
hi gh stresses in this region will be relieved by faulting, until the failure criterion is no longer 
satisfied. This region may thus be the most pervasively faulted and fractured zone, and thus 
perhaps a favored location for magma to collect. 

The issue of compressional versus extensional horizontal stresses in the volcanic edifice is 
critically important when considering volcano growth. In order for the volcano to grow, there 
must be a path for magma to reach the surface (or near the surface in the case of growth by 
intrusion). Large horizontal compressive stresses within the edifice may cut off summit eruptions 
unless there is sufficient magma overpressure. Thus, the route from a central magma chamber to 
the summit may be cut off, and flank eruptions would be favored, as noted above. After sufficient 
flank eruption, however, loading on the flanks might then lead to relief of compressive stresses 
near the summit, re-opening the path to the summit caldera. When the effects of time-dependent 
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aureole deposits arc of greatest extent to the north and west, but also there is a significant aureole 
lobe to the southeast. The basal scarp has greatest relief on the northwestern and southeastern 
flanks. Models with a detachment predict that at the edge of the edifice 03 = age and is 
extensional. This can account for radial normal faulting but not for the scarps and aureole 
deposits, which in the above scenario require that 03 = O n . The addition of a regional northwest- 
southeast extensional stress, induced by the regional gradient in elevation, however, would 
provided the necessary orientation of stresses in the northwestern and southeastern quadrants of 
the volcano, where the aureole and scarp are most prominent. The aureole deposits to the north 
and west of Olympus Mons could then be a result of large landslides analogous to those off the 
Hawaiian chain but on a larger scale. The basal scarp would then consist of the coalesced head 
scarps of these landslides. The lobate features to the northwest of Arsia and Pavonis Montes 
would then be similar in origin but on a smaller scale, without prominent basal scarp formation. 
This scenario also explains why the deposits are always more prominent on the northwestern side 
of these volcanoes, since this is the regional downslope direction. Note that the above scenario 
postulates a consistent stress field (of NW-SE extension) throughout the evolution of the volcano. 
This contrasts with the view of Francis and Wadge [1983], who suggest that proto-Olympus Mons 
was elongate to the northwest due to an early NE-SW extensional regime, which allowed 
preferential effusion along a NW-SE axis. 

Large landslides and slumps, triggered by extension on the flanks, may be a second important 
mode of volcano growth, in that mass movement clears the way for more volcanic material to be 
emplaced behind it [Lipman et al . , 1985]. The model of Lipman and others emphasizes the role of 
dike intrusion at the rift zones, with the subvolcano detachment passively taking up the resulting 
strain. Our results suggest that without such a detachment, the rift zones may not be able to form, 
due to a predicted environment of horizontal compression. As a welded volcano grows, an ever 
increasing horizontal compressive stress near the base of the edifice would eventually choke off 
magma supply to the summit. Thus, the Hawaiian rift zones may owe their existence to the 
presence of the detachment, and volcanoes lacking such a detachment may exhibit important 
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structural and evolutionary differences front Hawaiian volcanoes. For example, Nakamura [1980] 
has suggested Out the lack of rift rones in the Galapagos shields is anributable to a paucty of 

oceanic sediments on the young ocean floor beneath the volcanoes. 

Linear Rift Zones 

Unear symmetries observed in all the Tharsis Montes volcanoes may be indirect evidence for 
failure in die sub-volcano region. The second stage in to Cmnpler and Aubele [1978] sequence 
of volcano evolution involves die development of a linear rift bisecting die volcanoes, with 
eruptions emanating from this rift. Radially oriented normal faulting under a volcano nught be 
thought at first to tend to divide the volcano into radial sections, much like pie slices. The 
existence of the linear rifts on the Tharsis Montes, however, suggests that this mode of failure may 
be concentrated into one large linear feature that bUects tile volcano. This feature may owe its 
orientation to regional stress fields or to pre-existing zones of weakness. Our caiculations predict 
radially oriented normal faulting in a broad zone beneath a volcano. For tile case with T e = 20 km 
(figure 1 1 ), tins none extends beyond the volcano radius and reaches quite close to the surface. 

Such a zone, if concentrated into a linear feature and if continued to the surface by sustained 
faulting, may give rise to the linear rifts and bilateral symmetty observed on the Tharsis Montes. 
Stresses due to topography of Tharsis might be expected to add acomponent of northwest- 
southeast extension (the direction of greatest slope). This would concentrate the radial faulting mto 
a rift zone striking northeast-southwest, as is observed on tite Tharsis Montes. The rift zones 
would form orthogonal to the regional topographic gradient, as is observed. A problem with this 
mode, is that the rifting has to propagate upward somehow to the surface of tile lower volcano 
flanks The prediction of radially oriented normal faulting on the lower flanks for the model with a 
sub-volcanic detachment (Figure 20) suggests another mechanism for forming the observed rifts. 
This explanation has me advantage that rift formation would occur early, almost immediately after 
volcano creation, a timing more consistent with me chronology of Cnmpler and. Aubele [1978], If 
the detachment model is applicable to Mars, it may help explain the existence of the nfts. 
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The evolution of Pacific hotspot volcanoes may provide some insight into the formation and 
influence of a bisecting rift zone. For instance, Holcomb [1985] has suggested that the island of 
Molokai in the Hawaiian chain resembles a volcano cut almost exactly in half, with one part 
transported downslope in one or more massive submarine slides, and the other part remaining to 
form the present island. Bathymetry of the ocean floor north of Molokai [Moore, 1964] supports 
this conclusion. We have earlier conjectured that a similar but less extreme mechanism could 
explain the presence of graben on the upper slopes of Arsia Mons and Pavonis Mons, in that as the 
volcano halves slide apart from each other, an extensional environment is created near the summit. 
To investigate this possibility quantitatively, a fully three dimensional model, or an axisymmetric 

model with non-axisymmetric loading, is required. 

Detachment Tectonics 

The flank tectonic characteristics of the Tharsis volcanoes, combined with the observations of 
and modeling results for terrestrial volcanoes with detachments, strongly suggests that detachments 
analogous to those of Hawaii-type volcanoes [Lipman et al, 1985] exist on Mars. Such a 
detachment cannot arise from the presence of a thick oceanic sediment layer as for Hawaii. It has 
nonetheless been demonstrated that liquid water existed in moderate quantities at the Martian 
surface earlier in the history of the planet [e.g., Carr , 1987]. While liquid water is now unstable 
at the surface, large quantities of water are believed to be buried as ground-ice. Tanaka [1985] 
proposed that gravity sliding and spreading could account for the Olympus Mons aureole (and 
basal scarp) providing that a basal detachment zone of material with low shear strength was 
present. He calculated that a layer of material with about 10% interstitial or interbedded ice could 
provide the necessary detachment. The pre-volcanic upper crustal layer may have consisted 
dominantly of an impact breccia regolith produced by repeated impacts. Such a heavily fractured 
material is likely to be porous and thus may contain interstitial ice. Even in the absence of an 
excess ice content, the mechanical weakness of such a layer could contribute to the formation of a 
detachment between a strong lower crust and a volcanic edifice. Tanaka [1985] states that such a 
structure is unique to Olympus Mons. However, Zimbleman and Edgett [1991] have mapped 
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lobale deposits on the northwestern flanks of all three Th arsis Montes that may be smaller-scale 
analogues to the Olympus Mons aureoles, so some form of detachment beneath all four constructs 
appears to be indicated. 

Plate Thickness and Volcano Size 

The ratio of elastic plate thickness to volcano diameter is observed to have significant effects 
on the evolution of the subsidence, stress, and regions of predicted failure in our computations. 

We note several trends as plate thickness increases relative to volcano diameter. Study of the 
patterns of stress orientations and predicted failure in Figures 9, 1 1, and 13, shows that the outer 
region of failure, at a given time, increases in horizontal (and vertical) extent, and a greater 
proportion of the inner part of this region is predicted to fail at later stages by strike-slip faulting, 
with increasing plate thickness (Table 3). The deviatoric stress plots of Figures 10, 12, and 14 
show that the average depth and volume of the region of low deviatoric stress beneath the volcano 
also increase with increasing plate thickness. The magnitude of deviatoric stress within the 
volcanic edifice, and in particular the magnitude of maximum horizontal compressive stress, 
decreases with increasing plate thickness. As noted above, a zone of large horizontal compressive 
stress may affect the propagation of magma from its mantle sources to the volcanic edifice. Elastic 
plate thickness may thus play a role in determining the location of magma storage areas. For the 
largest plate thickness considered here, stresses in the edifice never exceed the failure limit (Figure 
13). 

The relative timing of events in our models is generally consistent with the chronology set 
forth by Crumpler andAubele [1978] for the Tharsis Montes. After construction of the main 
shield, rift formation and flank eruptions along this rift occur in the second stage of evolution. 
Figure 9b shows that, after 10 Maxwell times, the principal model result is the beginning of 
rotation of the stress axes in the edifice (which is expected to favor flank eruptions). After a few 
tens of Maxwell times the failure zone beneath the volcano (which could favor formation of the rift) 
has begun to grow. The third stage involves summit subsidence and formation of concentric 
fractures and graben. After 100 Maxwell times, our computations show increased subsidence as 


well as the growth of the region of normal faulting outward of the volcano (Figure 9c-g). The 
fourth stage of Crumpler and Aubele [1978) involves further eruption along the rift, flooding the 
«.,mmit and flanks. Late stage flank eruptions are consistent with the stress fields m our 
calculations. Summit eruptions, however, would not be favored by the computed stress fields. It 
is likely that the stress field is influenced by the effects of previous stages of the evolution. Flank 
eruptions result in an added flank load, not modelled here, to which the lithosphere will respond 
with additional flexure. Subsidence and additional flexure near the flanks of the volcano may at 
least partially relieve the compressions! stresses near the summit, thus enhancing the possibility of 
additional summit eruptions. By this scenario, the preferred site for eruptions could alternate 
between summit and flank, on time scales on the order of 100-1000 T M . Also, the effect of 
sections of the volcano sliding away from the bisecting rift may result in an extensional 
environment near the summit, once again favoring summit eruptions. 

Future Directions 

Because of the limitations of these calculations in accurately modeling the behavior of large 
volcanoes, we hope to be able to relax some of these limitation^ in future models. In our models, 
assume that the crust and mantle are perfectly viscoelastic; within an element, no provision is made 
for stresses that exceed reasonable failure limits (such as in all the shaded elements in Figures 9, 

11, and 13). The inclusion of localized yielding in a volcano model can be presumed to have 

important effects on the geometry and magnitudes of the calculated stress and deformation fields. 
A finite element code that incorporates plastic yielding will allow us to account more realisucally 
for the effects that failure will have on predicted fault locations and styles. Further, the effect of 
thermal stresses has been ignored. Thermal stresses can play an important role in the evolution of 
the volcano, both though stresses accumulated by cooling after emplacement of each layer, and 
through the time-dependent thermal stresses in regions reheated by the later passage or storage of 


magma. 
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CONCLUSIONS 

„ to paper we have used finite element calculations to sunulate the defonnational and stress 
response of the lithosphere to a volcano load. The results of these calculations are compared with 
observed tectonic features and infened eruption characteristics of large Mariian volcanoes. Some 
of rite zones of faUute and faulting predicted in our models have direct analogues in observed 
features; others may be indirectly related to observed structures. Three regions of failure with 
characteristic stress regimes are identified in our calculations. Flexure! stresses in tire lithosphere 
immediately surrounding the volcano can result in circumferential graben. HexuraUy reduced 
failure inawide zone beneath the center ofalarge volcano may play a major role re the 
development and modification of tire edifice. Given that magma ascending through tile lithosphere 
propagates along fractures oriented petpendicular to tire direction of leas, compressive stress, a 
rotation in tire principal stress orientations during lithospheric flexure in tire models tends to favor 
eruption sites shifting from the volcano summit to the flanks as flexure proceeds. We conclude 
*a, time-dependent lithospheric flexure is important in determining to location and style of 
faulting within and surrounding large volcanoes as well as regulating to timing and location of 

volcanic eruptions. 

Horizontal compressive stresses caused by flexure and subsidence may lead to ttnus, 
faulting on to flanks of large volcanoes. Incremental volcano growtii. however, results in a less 
extreme environment of horizontal compression on to upper flanks a, late stages of development. 
Lure load increments on to flanks are not stressed from eariier flexure, and horizontal stresses are 
no longer great enough to cause flank thrusting. The most highly stressed area of to edifice is 
near to top of to earliest load increment. Inclusion of a detachment between volcano and 
lithosphere results in an environment of horizontal extension in to pari of to volcano to, 
overthrusts to lithosphere. Without effects such as these, which reduce to extreme magnitude of 
horizontal compression in our mode, volcanoes, i, would he very difficult for a volcano to continue 

growth via summit eruptions. 



A mass movement origin has been suggested for the vast aureole deposits surrounding 
Olympus Mons [Harris, 1977; Lopes et al., 1980; Francis and Wadge, 1983] as well as smaller 
lobate deposits to the northwest of each of the Tharsis Montes [Zimbleman and Edgett, 1991]. 
Analogy with large slumps and slides off of the Hawaiian islands [ Moore et aL, 1989; Lipman et 
al., 1985] suggests that a detachment between volcano and lithosphere may play a role in the 
evolution of volc an oes on Mars as well as Earth. Models that include such a detachment, in 
contrast with models with welded volcanoes, predict a state of horizontal extension in the edifice 
(and eventually radial normal faulting on the lower flanks). When the effects of non-axisymmetric 
stresses due to regional topography are considered, predicted patterns of faulting match well 
observed structures on the Tharsis Montes. 
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FIGURE CAPTIONS 

Figure 1. Viking Orbiter view of Olympus Mons. The multiple-pit caldera and basal scarp (at 
bottom) are the main tectonic features visible. No graben are seen on the flanks of the volcano 
or on the volc ani c plains immediately surrounding the structure. Frame V0641A52, width of 

image is 440 km. 

Figure 2. (a) Viking Orbiter view of western Ascraeus Mons. The complex caldera is seen at the 
summit, and graben with a predominantly circumferential orientation are seen along the 
western and southwestern margins. Graben are located on Ascraeus Mons itself and in 
surrounding flows (unit as, Figure 3); some graben grade into collapse pits or sinuous 
channels formed by magma withdrawal and downslope flow. Frame V0643A78; width of 
image is 420 km. (b) Viking Orbiter mosaic of the western flank of Arsia Mons [Eliason et al. 

1991]. Note that graben extend all the way up the flank to the edge of the summit caldera, (c) 

Viking Orbiter mosaic of the eastern flank of Arsia Mons [Eliason et al. 1991]. Note the lack 
of graben and the radial texture of flows on the upper flank. Graben are evident near the 

boundary between the edifice and the surrounding plain. 

Figure 3. Geologic map of Ascraeus Mons and surroundings, simplified from Scott et al. [1981b] 
by Comer et al. [1985]. Volcanic units shown include relatively young Ascraeus Mons flows 
(as), intermediate-age Tharsis Montes flows (tm), and volcanic material undivided by age (vu). 
Also shown as a distinct unit is slide material (s), interpreted by Scott and coworkers as 
landslides and debris flows. Dashed lines show approximate elevation contours, m kilometers, 
relative to a fourth-degree, fourth-order equipotential [Wu, 1978]. The summit caldera 
complex is indicated by inward hatched lines. Extensional faults and graben are shown as 

heavy lines. 

Figure 4. Geologic map of Pavonis Mons and surroundings, simplified from Scott et al. [1981a, 
b, c] and Scott and Tanaka [1981] by Comer et al. [1985]. Units shown, in addition to those 
described in Figure 3, are relatively young volcanic flows from Pavonis Mons (pm) and Arsia 
Mons (am). Other information follows the format of Figure 3. Circumferential graben reach 
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quite far up the slopes, and an approximate bilateral symmetry can be seen about an axis 
trending approximately NNE-SSW. 

Figure 5. Geologic map of Arsia Mons and surroundings, simplified form Scott et al. [198 lc] by 
Comer et al [1985]. See Figures 3 and 4 for further explanation of symbols. Circumferential 

graben extend almost up to the summit caldera complex. 

Figure 6. (a) The finite element grid used for the calculations shown in later figures. The volcano 
is the small triangular area in the upper left-hand comer, (b) Plot of element types for models 
with incremental volcano loading, implemented by means of time-dependent element 
properties. Elements marked with a "3" are "on" from the start of the calculation; those with a 
"4" are switched "on" after 1000 Maxwell times, and those with a "5" are switched "on" after 

2000 Maxwell times. 

Figure 7. Comparisons of the TECTON solution and analytic thick-plate flexure [Comer, 1983] 
for vertical displacements resulting from an axisymmetric load simulating a conical volcano. 

For both cases, T e = tc = 40 km. (a) surface elevation relative to that at r=r max after 
essentially complete viscoelastic relaxation in the asthenosphere; (b) difference between the two 

solutions. 

Figure 8. Key for stress symbols in models. These symbols denote the principal stress directions 
under axisymmetry . The hourglass shapes are oriented along the direction of greatest 
compressive stress, the bars along the direction of least compressive stress. Magma is 
expected to ascend locally along the plane orthogonal to the least compressive stress, that is, 

perpendicular to the bar. 

Figure 9. Principal stress orientations and predicted failure in the vicinity of an incrementally 
grown volcano load on a 40-km-thick elastic lithosphere. See Figure 8 for the meaning of the 
stress orientation symbols. Elements in which stresses satisfy the Mohr-Coulomb failure 
criterion are shaded. Given the orientation of the stress symbol in a shaded element, the style 
and orientation of faulting can be inferred, (a) The elastic solution, (b-d) 20, 100, and 1000 
Maxwell times after the first load increment is added; (e-f) 20 and 1000 MaxweU times after the 
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second load increment is added; and (g) 1000 Maxwell times after the third load increment is 
added. 

Figure 10. Contours of maximum deviatoric stress in the vicinity ofthe volcano, forthe model of 

Figure 9. Contour interval is 50 MPa. (a) The elastic solution, (b-d) 1000 Maxwell times 
after the first, second, and third load increments are added, respectively. 

Figure 1 1 . Principal stress orientations and predicted failure for an incrementally grown volcano 
load on a lithosphere with T e = 20 km. (a) 1000 Maxwell times after the first load increment is 

added; (b) 1000 Maxwell times after the third load increment is added. 

Figure 12. Contours of maximum deviatoric stress for the model of Figure 11- (a) 1000 Maxwell 
times after the first load increment is added; (b) 1000 Maxwell times after the third load 

increment is added. 

Figure 13. Principal stress orientations and predicted failure for an incrementally grown volcano 
load on a lithosphere with T e = 60 km, 1000 Maxwell times after the third load increment is 

added. 

Figure 14. Contours of maximum deviatoric stress for the model of Figure 13, 1000 Maxwell 
times after the third load increment is added. 

Figure 15. Principal stress orientations and predicted failure for an incrementally grown volcano 
load on a crust 40 km thick and on an elastic lithosphere 60 km thick. 1000 Maxwell times after 
the third load increment is added. 

Figure 16. Contours of maximum deviatoric stress for the model of Figure 15, 1000 Maxwell 
times after the third load increment is added. 

Figure 17. Principal stress orientations and predicted failure for an instantaneous volcano load on 
a lithosphere with T c = 40 km, after 1000 Maxwell times. 

Figure 18. Contours of maximum deviatoric stress for the model of Figure 17, after 1000 

Maxwell times. 

Figure 19. Principal stress orientations and predicted failure for a model with a growing volcano 
load, a basal detachment between lithosphere and volcano, and crest and elastic lithosphere 40 
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km thick, (a-b) 10 and 1000 Maxwell times after the first load increment is added; (c) 1000 

Maxwell times after the third load increment is added. 

Figure 20. Contours of maximum deviatoric stress for the model of Figure 19. (a) 1000 MaxweU 
times after the first load increment is added; (b) 1000 MaxweU times after the third load 

increment is added. 

Figure 21. Principal stress orientations and predicted failure for a growing volcano load with 
isotropic prestress in the lithosphere. Crust and elastic lithosphere are 40 km thick, (a-b) 40 
and 1000 Maxwell times after the first load increment is added; (c) 1000 Maxwell times after 

the third load increment is added. 

Figure 22. Mohr diagrams for typical elements in regions of predicted failure. Diagonal line is the 
Mohr-Coulomb failure envelope for basalt [Handin, 196 6] with an angle of internal friction <j> = 
49° and cohesion c = 38 MPa. (a) Mohr diagram for an element at the top of the lithosphere 
(stresses calculated at element center r = 265 km, z = -2.5 km). In order of increasing radius, 
the circles represent the stress state at t = 0 (elastic solution), 20, 100, 1000, 2000, and 3000 
x M . When the failure envelope is first exceeded (r = 70 x m ), 03 = a n and <3\ = <?zz- At the 
final time shown, o 3 = o n and a, = a 00 . See Figure 9 for the principal stresses used to 


calculate circles. 


TABLE L Dimensions of Tharsis Volcanoes 


Volcano 

Diameter, 

km 

Relief, 

km 

Olympus Mons 

600* 

24 a 

Ascraeus Mons 

400 b 

18 b 

Pavonis Mons 

320 b 

14 b 

Arsia Mons 

420 b 

19° 


a Wu et al. (1984). 
b Blasius and Cutts (1976). 
c Blasius and Cutts (1981). 
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TABLE 2. Adopted Parameter Values 
Parameter Crust Mantle 

E, Pa lxlO 11 3xlO n 

p, kg/m 3 3000 3500 

v 0.25 0.25 


Lithosphere Asthenosphere 
tl, Pa s lxlO 27 lxlO 21 


? 
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TABLE 3. Characteristics of Predicted Near-surface Failure for Each Model. 


Model 1 / 
increment 

Inner radius 2 , 
km 

Outer radius 2 ’ 
strike-slip 
faulting, km 

Outer radius 2 * 
normal faulting, 
km 

First observed 
failure (time, 
radius, type) 

GW40 1 

140 

230 

380 

20 tm 

2 

150 3 

270 

440 

130-170 km 

3 

170 3 

300 

490 

mixed 

GW20 1 

120 

170 

280 

30 tm 

2 

120 3 

190 

320 

120-140 km 

3 

140 3 

220 

360 i 

radial normal 

GW60 1 

180 

270 

410 

40 x M 

2 

180 

330 

520 

170-220 km 

3 

190 3 

360 

590 

mixed 

GW40/60 1 

200 

300 

430 

40 tm 

2 

200 

350 

550 

200-220 km 

3 

210 

390 

630 

radial normal 

IW 40 1 

190 3 

290 

490 

10 TM 





200-250 km 





radial normal 

GD40 1 

120 

240 

380- 

20 tm 

2 

. 130 3 

280 

440 

120-150 km 

3 

150 3 

310 

500 

mixed 

GWPS40 1 

150 

290 

340 

40 tm 

2 

160 

310 

410 

140-190 km 

3 

170 3 

320 

470 

mixed 

1 Model key: 



G = growing volcano, 

W = welded volcano, 

1 = instantaneously emplaced volcano, 

D = detached volcano, 

PS = isotropic pre-stress applied to lithosphere, integer = elastic lithosphere 
thickness. 

2 Radius values given at t = 1000 tm after application of given load increment. 

3 Value less than the maximum radius of the volcano at a given load increment 
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Arsia 2.2 13 40km time= 1.66700E4-11 seconds 

stress orientation/failure plot full output timestep # 4 

' cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Arsia 2.2 13 40km time= 8.33300E+12 seconds 

stress orientation/failure plot full output times tep # 12 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Area 2.2 13 40km time = 1.66700BH3 seconds 

stress orientation/failure plot full output timestep# 22 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Figure 9f 






































































































































































Arsia 2.2 13 40km 

stress orientation/failure plot ftjlounmttimestep# 

cohesion = 3.50000Bf07 pa; phi = 49.0000degrees 
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Figure 11a 





























































Aim 2.2 13 20km *”e = 2.50000Bt-13 seconds 

stress orientation/failure plot full output timestep # 32 

cohesion- 3.5OOOOB*07 pa; phi = 49.0000degrees 
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Arsia 2.2 13 60km time = 8.33300E+12 seconds 

stress orientation/failure plot full output times tep # 12 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Figure 14b 






Arsia 2.2 13 crust=40km lith=60km time= 2.50000E+13 seconds 

stress orientation/failure plot ^ output times tep # 32 

cohesion* 3.50000E+07 pa; phi = 49.0000degrees 
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Figure 15 
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Figure 17 





Arsia 2.2 13 slip 40km time = 8.33300E+10 seconds 

stress orientation/failure plot fall output times tep # 2 

cohesion - 3.50000B+07 pa; phi = 49.0000degrees 
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Arsia 2.2 13 slip 40km **“ = S.33300BH2 seconds 

stress orientation/failure plot full output timestep# 12 

cohesion = 3.50000B+07 pa; phi — 49.0000degrees 
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Figure 19b 
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Arsia 2.2 13 40km psl time = 3.33300E+1 1 seconds 

stress orientation/failure plot full output times tep # 6 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Figure 21a 







































































































































































































































Arsia 2.2 13 40km psl time = 8.33300E+12 seconds 

stress orientation/failure plot full output times tep # 12 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Figure 21b 
































































































































































































































Arsia 2.2 13 40km psl time = 2.50000E+13 seconds 

stress orientation/failure plot faU output times tep # 32 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Features on Venus Generated by Plate Boundary 
Processes 

Dan McKenzie (Bullard Laboratories, Madingley Road, Cambridge, 
CB3 OEZ, U.K.; ph. 44-223-337177; fax. 44-223-60779; e-mail 
dpml @phx.cam.ac.uk); Peter G. Ford (Center for Space 
Research, MU, Cambridge, MA 02139); Catherine Johnson and 
David Sandwell (Scripps Institution of Oceanography, La Jolla, 
CA 92093); Barry Parsons (Department of Earth Sciences, 
Oxford University, Oxford 0X1 3PR, U.K.); Stephen Saunders 
(Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, 

CA 91 109); Sean C. Solomon (Dept, of Earth, Atmospheric and 
Planetary Sciences, MU, Cambridge, MA 02139) 

SAR images of Venus, taken with a radar whose wavelength is 12.6 
cm, are compared with GLORIA images of active plate boundaries, 
obtained with a sound source whose wavelength is 23 cm. Features 
similar to transform faults and to abyssal hills on slow and fast 
spreading ridges can be recognized within the Artemis region of 
Venus, but are not clearly visible elsewhere. The composition of the 
basalts measured by the Venera 13 and 14 and the Vega 2 spacecraft 
corresponds to that expected from adiabatic decompression, like that 
which occurs beneath spreading ridges on the Earth. Structures that 
resemble trenches are widespread on Venus, and show the same 
curvature and asymmetry as they do on Earth. These observations 
suggest that the same simple geophysical models that have been so 
successfully used to understand the tectonics of the Earth can also be 
applied to Venus. 


THE GABBRO - ECLOGITE PHASE TRANSITION AND THE 
ELEVATION OF MOUNTAIN BELTS ON VENUS; Noriyuki Namiki and Scan C. 
Solomon, Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of 
Technology, Cambridge, MA 02139. 

Introduction. The linear mountain belts of Ishtar Terra on Venus are notable for their 

topographic relief and slope and for the intensity of surface deformation [1,2]. The mountains 
surround the highland plain Lakshmi Planum. Volcanism is rue to absent in Maxwell, Freyja, and 
A Win Montes, but a number of magmatic features are evident in Danu Montes [2,3], the mountain 
range least elevated above Lakshmi Planum. Whether western Ishtar Terra is a site of mantle 
upwelling and consequent hot spot volcanism [4-6] or of mantle down welling and consequent, 
convergence of lithospheric blocks [7,8] is currently a matter of debate. However, the mountains 
are generally regarded as products of large-scale compression of the crust and lithosphere I2»yj. 

Among the four mountain belts surrounding Lakshmi Planum, Maxwell Montes is the mghest 
and stands up to 1 1 km above the mean planetary radius and 7 km above Lakshmi Planum. The 
bulk composition and radioactive heat production of the crust on Venus, where measured, are 
similar to those of terrestrial tholeiitic basalt [10]. Because the thickness of the low-density OTist 
may be limited by the gabbro - garnet granulite - eclogite phase transitions (Fig. 1), the 7-1 1 km 
maximum elevation of Maxwell Montes is difficult to understand except in the unlikely situation 
that the crust contains a large volume of magma [1 1]. A possible explanation is that the base of the 
crust is not in phase equilibrium. It has been suggested that under completely dry conditions, the 
gabbro - eclogite phase transition takes place by solid state diffusion and may require a geologically 
significant time to run to completion [12]. Solid state diffusion is a strongly temperature-dependent 
process. In this paper we solve the thermal evolution of the mountain belt to. attempt to constrain 
the depth of the gabbro - eclogite transition and thus to assess this hypothesis quantitatively. 

Thermal Model. The one-dimensional heat equation is solved numerically by a finite difference 
approximation. The deformation of the horizontally shortening crustal and mantle portions, of the 
thermal boundary layer is assumed to occur by pure shear, and therefore the vertical velocity is 
given by the product of the horizontal strain rate V and depth z. The thermal diffusivity is assumed 
to be l.OxlO^mV 1 in both crust and mantle. Crustal heat production is assumed to equal 1.4x10* 
13 k s* 1 . The initial temperature profile is determined by the assumption of steady-state conditions 
with zero velocity. Temperature at the surface and the bottom of the thermal boundary layer are 
fixed at 750 K and at a value T u taken as free parameter. .. ... 

The phase dia g ram is assumed to be that of tholeiinc basalt [13], and the densities of gabbro 
and eclogite are taken to be 2900 and 3500 kg nr 3 . The density of gamet-granulite is assumed to 
increase linearly from that of gabbro to that of eclogite as pressure increases at a given temperature. 
The density of the mantle is assumed to be 3400 kg m* 3 . The micromechanism of the gabbro - 
eclogite transition is not well understood. In this study we assume that the volume diffusion of 
cations is the most likely rate-limiting process of the transformation, which involves chemical as 
well as phase changes. The volume fraction of reacted component, y, is given by, 

\jn|/ = D/r^ 

where r is the grain radius and D is the diffusion coefficient [12). Since the slowest-moving cation 
limi ts the reaction rate and AP+ is likely to be this cation, we adopt as a minimum value for D the 
lower end of the range of estimates of the diffusion rate of Al 3 * in orthopyroxene [14]: 

D = Dai, Opx = l.lxl0' 5 exp (-400 kJ/RT ) 

where R is the gas constant and T is the absolute temperature. "Hie diffusion rate, however, is 
experimentally uncertain because AP+ diffusion is extremely sluggish, particularly at low 
temperature. In order to bound D from above, we use the diffusion rate of Fe in garnet (L>Fe.Gt - 
6.4x1 0* 8 exp [-( 270 kJ + 5.7xl0* 5 P )/RT]), where P is in Pa. For each parcel of shortening 
lithosphere, \jr is obtained by integration over time. The density at a given depth is determined from 
the volume fractions of unreacted and reacted components. , 

Numerical Results. Temperatures in the thickening crust and mantle are calculated for rates of 
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horizontal convergence of lO* 1 ^ (Fig. la) and lO* 1 ^ s* 1 (Fig. lb). For all models discussed here 
(Table 1), thicknesses of crust and thermal boundary layer are assumed to be initially 20 and 50 
km, respectively, and to increase to 100 and 250 km, respectively. Temperature profiles for the 
strain rate of 10* 15 s* 1 are vertically stretched as the crust and lithospheric mantle are thickened 
(Fig. la). Temperatures do not increase significantly from initial values because heat is mainly 
transferred by advection and the contribution of crustal heat production is minor. Hence gabbro 
remains metastable for 50 My or more, and the elevation of shortened lithosphere can increase as 
much as 12 km above the surrounding plains in that time interval (Fig. 2a). The phase transition 
proceeds, i.e., elevation is limited, only if grains are small and diffusion is fast (Model 3). 

For Y=10 _1 ® s* 1 , crustal heat production dominates advective heat transfer after the crust 
becomes as thick as 60-80 km (Fig. lb). The resulting increase in temperature hastens die phase 
transition. The slower strain rate also lengthens the formation time of the mountains relative to the 
characteristic reaction time, which depends on r, D, and temperature at the base of the crust For 
larger grains (r=10 mm), the elevation reaches 1 1 km or more if the initial temperature at the base 
of the crust T c , is 1 150 K (Model 6 in Fig. 2b). For the same value of T c the elevation is at most 6 
km for grains of 1 mm radius (Model 7). This result constitutes an upper bound on T c for small 
grain size. If D=DFe,Gt is assumed, that upper bound is lowered to 1050 K (Model 8). 

Discussion. Because at long wavelengths the topography of western Ishtar Terra is correlated 
with the gravity field, dynamical support of the broad 4-km elevation of the region is likely [e.g., 
15). Therefore the 7-km elevation of Maxwell Montes above the adjacent plateau is a more 
meaningful constraint on maximum relief than the 1 1-km elevation above mean planetary radius. 
The results are insensitive to the assumed initial thickness of crust but are sensitive to the density 
difference between crust and mantle. If densities of 3000 and 3300 kg m* 3 for crust and mantle are 
flgqimftH, elevations are -40% lower than the values presented here. Such changes constrain the 
mo del s of thermal structure and phase transition depth more severely at low strain rate than the 
density values adopted above. 

Two diffusion rates have been assumed in this study so as to represent a wide range of 
diffusion data. We should note, however, that under wet conditions, i.e., in the presence of either 
water [16] or melt [17], grain-boundary diffusion becomes much more efficient than volume 
diffusion. This is potentially noteworthy for understanding the contrast between Maxwell and 
Danu Montes. Despite the fact that Danu Montes display compressional deformation as extensive 
as the other mountain belts, the maximum elevation is as little as 0 km above the bounding plateau. 
Such comparatively modest elevation may be related to the presence of magmatic features within 
Danu Montes, if elevation is limited by an enhanced diffusion rate because of the melt at grain 
boundaries. Assessing the cause of higher temperatures beneath Danu Montes requires more 
detailed thermal models than the simple one-dimensional model considered here. 

Conclusions. Taking into account the temperature-dependent reaction rate of the gabbro - 
eclogite phase transition, horizontal strain rates of 10* 15 and 10* 16 result in significant differences 
in the of maximum elevation of mountains, not only because of the difference in the formation time 
for relief, but also because of the difference in the thermal regime from advection-dominanted to 
crustal-heat-production dominated. For V^IO* 15 s* 1 , the observed maximum elevation of mountain 
belts can be explained as a consequence of disequilibrium phase boundary depth for a wide range 
of physical parameters, although a comparatively young age for Maxwell Montes (50 My) is 
implied. For the lesser horizontal strain rate of 10* 16 s* 1 , only limited parameter values for thermal 
models are allowed. 

References. [1] V.L. Barsukov et al., JGR, 91, D378, 1986; [2] S.C. Solomon et al., 
Science, 252, 297, 1991; [3] J.W. Head et al.. Science, 252, 276, 1991; [4] A. A. Pronin, 
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Abstract 

Maxwell Montes, standing up to 7 km above the adjacent highland plateaus, constitute the 
hi ghftgf mountain belt on Venus. Because the thickness of the crust is likely to be limited by the 
gabbro - garnet granulite - eclogite phase transitions, this relief is difficult to reconcile with the 
assumption of thermodynamic equilibrium and a standard Airy isostatic model. We explore the 
hypothesis that the crust-mantle boundary is not in phase equilibrium, but rather is rate-limited by 
the temperature-dependent diffusion of the slowest ionic species. Under the simplifying 
assumption that the mountains formed by uniform horizontal shortening of the crust and 
lithospheric mande at a constant rate, we solve the one-dimensional thermal evolution problem. 

The time-dependent density structure and surface elevation are calculated by assuming a 
temperature-dependent reaction rate and local Airy isostatic compensation. For a rate of horizontal 
strain of 10*^ s'^ or greater, the rise in temperature at the base of the crust during mountain 
formation is modest to negligible, the deepening lower crust is metastable, and surface elevation 
increases as the crust is thickened. For strain rates less than 10'*^ s'^, in contrast, crustal 
temperature increases with time because of internal heat production, and the lower crust is more 
readily transformed to the dense eclogite assemblage. For such models a maximum elevation is 
reached during crustal shortening. While this maximum relief is 7 km or more for some models, a 
smaller density contrast between crust and mantle than assumed here (500 kg nr 3 ) and 
incorporation of horizontal heat transport would lessen this value. We therefore favor formation of 
the mountain belt at a strain rate at least of order 10- 15 r 1 . By this reasoning. Maxwell Montes 
must be comparatively young, of order 50 My. 
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INTRODUCTION 

The linear mountain belts of western Ishtar Terra on Venus are notable for their topographic 
relief and slope and for the intensity of surface deformation [Barsukov et al., 1986; Solomon et al., 
1991, 1992]. These four mountain belts, including Akna, Freyja, Maxwell, and Danu Montes 
(Figure 1), are generally regarded as products of large-scale compression and horizontal shortening 
of the crust and lithosphere [Crumpler et al., 1986; Solomon etal., 1991, 1992; Kaula et al., 

1992]. Among these mountain belts. Maxwell Montes is die highest and stands up to 1 1 km above 
mean planetary radius and 7 km above the highland plain Lakshmi Planum (Figure 1). Volcanism 
is rare to absent in Maxwell, Freyja, and Akna Montes, but a number of magmatic features are 
evident in Danu Montes [Solomon et al., 1991; Head et al., 1991; Kaula et al., 1992], the 

mountain range least elevated above Lakshmi Planum. 

The bulk composition and radioactive heat production of the crust on Venus, where measured, 
are similar to those of terrestrial basalts [Surkov et al., 1984, 1987]. Because basalt transforms to 
eclogite at high pressure and because eclogite is likely to be denser than mantle material, the 
thickness of the crust may be limited by the gabbro - garnet granulite - eclogite phase transitions 
[Anderson, 1981; Turcotte, 1989]. The phase transition depth depends on the temperature 
structure of the crust beneath the mountains. If the thermal gradient is low and thermodynamic 
equilibrium is assumed, the phase changes would take place at shallow depth, and simple isostatic 
models underpredict the relief [V order Bruegge and Head, 1991]. If the thermal gradient is high, 
in contrast, the base of the crust is predicted to undergo melting before a sufficiently deep root can 
form to support the mountains. Thus the 7-1 1 km maximum elevation of Maxwell Montes is 
difficult to understand under the assumptions of local isostasy and thermodynamic equilibrium 
except in the unlikely situation that the crest beneath the mountains contains a large volume of 
magma [Vorder Bruegge and Head, 1991]. 

A possible explanation for the great relief of Maxwell Montes is that the base of the crust is not 
in phase equilibrium. Because of the high surface temperature on Venus (750 K) and the very low 
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water abundance of the lower atmosphere [von Za/w er a/., 1983], the Venus crust is thought to 
contain negligible water [Kaula, 1990]. It has been suggested that under anhydrous conditions, 
the gabbro - eclogite phase transition takes place by solid-state diffusion and may require a 
geologically significant time to run to completion [Ahrens and Schubert , 1975]. Solid-state 
diffusion is a strongly temperature-dependent process, so that quantifying this suggestion for 
application to Venus requires the solution of a heat transport problem. In this paper we develop 
simple models for the thermal evolution of the crust beneath mountain belts on Venus in an attempt 
to constrain the time-dependent depth of the gabbro - eclogite transition and thus to assess this 
hypothesis. 

Because at long wavelengths the topography of western Ishtar Terra is correlated with the 
gravity field, at least partial dynamical support of the regional elevation of about 4 km is likely 
[e.g., Grimm and Phillips, 1991]. We therefore adopt the 7-km elevation of Maxwell Montes 
above the adjacent plateau as a conservative estimate of the relief attributable to local isostasy. 

THERMAL MODEL 


We assume that the crustal and mantle portions of the thermal boundary layer shorten 
horizontally and thicken vertically with time in the manner of pure shear. The thermal structure is 
then governed by the one-dimensional heat equation. 


ar + , z 9r = K ^ + ^ 


( 1 ) 


dt ' _ dz Bz 2 Cp 
where T is temperature, t is time, yis the horizontal strain rate, z is depth, k is the thermal 
diffusivity, A is the crustal heat production, and Cp is die specific heat. Under the assumption that 
deformation is by pure shear, y is constant and the vertical velocity is given by the product of y and 
z. We solve equation (1) numerically by an explicit finite difference approximation. 


Model Parameters 

The thermal diffusivity is assumed to be l.OxlO- 6 m 2 s 1 in both crust and mantle. The crustal 
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heat production is assumed to equal l.OxlO- 10 W kg* 1 , on the basis of the K, U, and Th 
concentrations determined by gamma ray spectrometry by Venera and Vega landers [Surkov et al., 
1987]. The specific heat is assumed to be 850 kJ kg* 1 K* 1 . The other physical properties of 
crustal and mantle materials are assumed to be those, respectively, of tholeiitic basalt and peridotite 
[Basaltic Volcanism Study Project, 1981]. Temperature at the surface and the bottom of the 
thermal boundary layer are fixed, respectively, at 750 K and at a value T w , a free parameter in 
numerical models. The initial temperature profile is determined by the assumption of steady-state 
conditions with zero strain rate. The increase in temperature with depth due to adiabatic 
compression is not considered here, because compressional heating is negligible compared with 
crustal heat production. 

Tj>i is interpreted as the temperature in an isothermal core of a thermal convection cell [Turcotte 
and Schubert, 1982] and is, thus, fixed at a constant For a given Tbi, an initial temperature T c at 
the base of the crust is determined. Since temperature at the base of the crust controls the elevation 
of the mountains, T c has more importance than T w . Zuber [1987] analyzed wavelengths of 
tectonic features on Venus and found that the crustal thickness is less than about 15 km if the 
average vertical thermal gradient in the crust dT/dz is 25 K km* 1 , and less than about 30 km if 
dT/dz is 10 K km* 1 . The absence of significant viscous relaxation of impact crater relief also 
limits the crustal thickness in plains regions to be less than 10 km for dT/dz = 20 K km* 1 and 20 
km for dT/dz = 10 K km* 1 [Grimm and Solomon, 1988]. These results give a maximum 
temperature difference across the crust of about 400 K. Therefore T w is given initially so that T c 
satisfies this upper bound. 

The phase diagram is assumed to be that of tholeiitic basalt [Ito and Kennedy, 1971], and the 
densities of gabbro and cclogite arc taken to be 2900 and 3500kg nr 3 . The density of garnet 
granulite is assumed to increase linearly from that of gabbro to that of cclogite as pressure increases 
at a given temperature. The density of the mantle is assumed to be 3400kg nr 3 . The elevations 
calculated by an assumption of local Airy isostatic compensation arc sensitive to the density 
difference between crust and mantle. We have adopted a relatively large value of 500 kg m*3. 


Because a large density difference results in a higher elevation for a given root thickness, this 
assumption permits a conservative lower bound on horizontal strain rate. We discuss later the 
consequence of a lesser density contrast 

Diffusion Rates 

The micromechanism of the gabbro - eclogite transition is not well understood. In this study 
we assume that volume diffusion of cations is the process most likely to limit the rate of the 
transformation, which involves chemical as well as phase changes. The volume fraction of reacted 
components, Y. is given by, 

YY = D/r 2 (2) 

where r is the typical grain radius and D is the diffusion coefficient [Ahrens and Schubert , 1975]. 
For each parcel of shortening thermal lithosphere, Y is obtained by integration over time. The 
density at a given depth is determined from the volume fractions of reacted and unreacted 
components. 

There are several reactions occurring among constituent minerals in the garnet granulite stability 
field [Ringwood, 1975]. The rate of each reaction is difficult to assess, because of a paucity of 
experimental da t a We focus on the diffusion rates of ionic species in pyroxene, garnet, and 
olivine, because those minerals characterize the phase assemblages and therefore play important 
roles in the phase transitions. Plagioclase is possibly another important mineral for this 
assessment, but no data for self diffusion in plagioclase are available. The data for volume 
diffusion in silicates show that silicon and oxygen diffuse much more slowly than divalent and 
trivalent cations. Thus the Si-O groups provide a static framework through which the cations 
diffuse [Freer, 1981]. We examine the diffusion rates of four major cations, Mg 2+ , Fe 2+ , Ca 2+ , 
and Al 3 * (Figure 2). The slowest diffusion rate is likely to limit the reaction rate. 

Diffusion rates of Fe 2+ and Mg 2+ in garnet and olivine [Chakraborty and Ganguly, 1991; 
Morioka and Nagasawa, 1991], of Mg 2+ in pyroxene [Cygan and Lasaga , 1985], and of Ca 2+ in 
olivine [Morioka and Nagasawa, 1991] have been experimentally determined (Figure 2). Of these 
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the diffusion rate of Mg 2+ in garnet (DMg,Gt) is dte lowest Of course. , ^ diffusions of aey 
unmeasured major species is slower than that of Mg 2 * in garnet, then the phase change rate will be 
limited by this slower cation. We therefore take the diffusion coefBoent t> Mg .G. for thts process as 

an upper bound on D in equation (2). 

Other diffusion rates, while no, welMetermined experiment, have been esdmared, either 
fiom measurement at a single temperance, or from observedcompositional gradients. Whde all 
such estimates include large uncertainties, we take die minimum value of these estimates as a lower 
hound on D in equation (2). This minimum is for diffusion of A1 in orthopyroxene (Figme 2), 

estimated from an analysis of compositional gradients in natural assemblages ^ Yhej”' 

,991] The estimated rate depends on assumptions of cooling rate and pressure. Whdenetther 

dam nor estimates are available for die diffusion coefficients of Al 3 * in olivine and 

diffusion in olivine and garnet is generally faster than in pyroxene [Freer, i981; Smithand Brown, 

1988], and thus is no, likely to be rate-limiting. The diffusion of Ca 24, in garnet has been estimated 

by numerical sinndation only at 1200*C [ChakrabortyandGangufy, 1991] and is 

AP» in orthopyroxene at that temperature (Figure 2). The mobilmes of Ca and 

be comparable in gamet and pyroxene, because neither garnet nor pyroxene in natural assem ges 

is consistently homogeneous in those elements [Smith and Barron, 1991]. These considerations 

reinforce the view that DAhOpx gives a reasonable lower bound on D in equation (2). The two 
bounds on D are 

D £ D Mt ,o. = 2.8 x 10 4 exp [- (270 kJ + 3.2 x Kh 6 P )/( RT )J 
D SDai. Op* = 1-1 x 10* 5 exp 1-400 kJ /( RT )] 

where Pis in Pa [Moriofci and Nagasawa, 199U Smith and Barron, 1991]. 

There is some possibility that die assumed lower bound on D is too large. Sooner « of. [ 98 
and Sautter and Harte [1990] reported a diffusion rare of Al 2 * in clinopyroxene less than the lower 
hound given above (Figure 2), while Freer e. of. [1988] experimentally obtained a maximum 
diffusion rare as great as die upper bound given above (Figure 2). Further, die diffusion rates in 
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plagioclase are generally less than in pyroxene [Smith and Brown, 1988], and interdiffusion of 
CaAl-NaSi in plagioclase is several orders of magnitude slower than die lower bound on D [Grove 
et al., 1984]. If D is as low as suggested by these last results, no phase transition is expected to 
occur during the formation of mountains on Venus. However, the reaction rate is unlikely to be 
much slower than as given by the lower bound on D, because if volume diffusion is too slow, the 

reaction rate will be limited by grain-boundary diffusion. 

Grain-boundary diffusion, in fact, has been proposed as a possible governing process m 
metamorphic reactions, because grain-boundary diffusion coefficients can be several orders of 
magnitude greater than those for volume diffusion [Joesten, 1991]. Whether the phase change 
proceeds by volume diffusion or grain-boundary diffusion will be governed by the process that is 
most efficient for mass transport Grain-boundary diffusion coefficients can be misleading, 
however, because mass transport by such a process is confined to a thin layer, i.e., the grain 
boundary. Unfortunately, grain-boundary diffusion data for the minerals involved in the gabbro - 
garnet granulite - eclogitc reactions are few and generally not well determined [Joesten, 1991]. 
Reliable experimental data for grain-boundary diffusion of oxygen in forstcrite, for instance, have 
been reported only at two temperatures [Condit et al., 1985]. These measurements suggest that 
grain-boundary diffusion for grain sizes of 1-10 mm is no more efficient for bulk transport than the 
lower bound on volume diffusion adopted above (Figure 2). 

NUMERICAL RESULTS 

Temperatures in the thickening crust and mantle are calculated for rates of horizontal 
convergence of 10-15 (Figure 3a) and 10*“ s-> (Figure 3b). For all models discussed here (Table 
1), thicknesses of crust and thermal boundary layer are assumed to be initially 20 and 50 km, 
respectively, and to increase to 100 and 250 km, respectively. Temperature profiles for the strain 
rate of 10- 15 s 1 are vertically stretched as the crust and mantle portion of thermal boundary layer 
are thickened (Figure 3a). Temperatures do not increase significantly from initial values because 
heat is transferred mainly by advection and the contribution of crustal heat production is minor. 
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Hence gabbro remains metasrabl. for 50 My or more, and An elevation of shonened lithosphere 
can increase as much as 12 km above the sunounding area of undeforaned lithosphere in that time 
interval (Model 1 in Figure 4a). Equation (2) shows that the reaction runs to completion in a 
shorter time for small grains than for larger grains. Also it is clear from equation (3) that a higher 
temperature increases the diffusion rate and promotes the reaction. For the slower diffusion ra 
(PMCt). however, a strain rate of 10-'5 s-' is sufficiently to to suppress these effects (Model 2). 
The phase transition proceeds. i.e„ elevation is limited, only if grains me small, diffusion is rapid. 

and temperature is high (Model 3). 

FUr a strain rate of 10-'« s ', crustal heat production dominates adveedve heat transfer after the 
crust becomes as thick as 60-80 km (Figure 3b). The resulting increase in temperature hastens the 
phase transition, and the slower strain rate lengthens the formation time of the mountains. For 
larger grains 0=10 mm) and a diffusion constant given by equation (3a) the elevation reaches 1 1 
km or more even if T e is as great as 1 120 K (Model 4 in Figure 4b). For the same value of T c the 
elevation is at most 7 km for grains of 1 mm radius (Model 5). If T c is greater than 1120 K, the 
model gives a maximum elevation lower than 7 km, that is, the model fails to explain the observed 
relief of Maxwell Montes. Therefore this result provides an upper bound on T c for small grain 
size. If D is given by equation (3b), that upper bound is lowered to 1030 K (Model 6). Thus the 
combination of the longer formation time due lo the lesser value off and the faster reaction due to 
higher temperature yields constraints on the range of parameters controlling the phase change if the 
elevation of the mountains is to reach 7 km or more above the adjacent plateaus. 

DISCUSSION 

The thermal model described above is simple, but significant constraints on the strain rate 
governing mountain building and the temperature at the base of the crust are nonetheless obtained. 
In this section, we firs, consider how the results would differ if key parameters were varied or the 
models were modified. We then discuss some implications of the numerical results given above. 

The thermal evolution results are insensitive to the assumed initial thickness of the crust At a 


V I'll' I """' 1 
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« of ID-' V, .he advection-dominared thermal profile >•*— to the toperato of 
each parcel of crus, changes hole as to parcel deepens. For a given "due of T c . «be thermal 
structure ar a given crustal thickness is not affected by initial crustal thickness, or, equivalently, 
initial thermal gradient Only to time required for thecrcstto reach to tohness depends on the 
initial thickness. A, a strain rate of 10" 16 s* 1 , in contrast, die thermal profile is close to steady swte, 
to is, the geothertn is neatly die same as to for zero strain rare and thus does not depend on the 
initial crustal thickness. The into cmstal thickness is sdU impexto. for the maximum rehef, 
however, because by assumption it is identical to the crustal thickness beneath surrcundutg areas, 
taken here to be undergoing negligible rates of horizontal strain. By the assumption of Any 
isostasy, at a given cmstal thickness beneath the mountains a difference of 10 km m the to ess 

of the crust beneath surrounding areas results in a 1.5 km difference in rehef. 

Also important is the density difference between crus, and mantle. If densities of 3000 an 
3300 kg nr 3 for crust and mantle are assumed, relative elevations are -40% lower than the values 
presented here. A, a strain rate of 10 45 s' 1 , gabbroic lower crust remains metastable during die 
formation of mountains in most cases, as shown by models 1 and X Therefore a smaller density 
contras, simply requires thicker crost, i.e„ greater duration of shoeing, to does no, otherwise 
constrain to model parameters. In contrast, a, a strain rare of 10-« s-'. maximum elevauons are 
limited in most cases by to gabbro - eclogire phase transition because of to high temperature a, 

the base of to crust, as shown by models 5 to 6. The lower relief resulting horn a smaller 

density contrast lessens to upper bound on T c . Thus such a lower density contrast constrains tire 
models of thermal structure and phase transition depth more severely at low strain rates than do to 

density values adopted above. 

A, a strain rate of 10 " 16 s' 3 , to thermal regime is dominated by crustal heat production and 
diffusion of heat. Horizontal hea, transfer within to shorrening thermal lithosphere therefore re 
probably not negligible. The maximum increase in, emperature at to totiom of to crus, due to 
horizontal hea, transfer can be estimated horn to solution for to insrentaneous hearing of a semt- 

infinite half-space [Turcotte and Schubert, 1982], 
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ST^Tcrfcj^) 


(4) 


where ST is the increase in temperature, AT is the temperature difference beneath die crust and an 
isothermal core which surrounds the shortening lithosphere, L is the horizontal scale of the 
mountains, and t is the formation time. In our problem, AT is given by Tj>i - T c , and L is -100 
km. Formation times of 50 My at a strain rate of 10*15 s*l and of 500 My at a strain rate of 1046 
s* 1 result in ST of 2xl0* 2 and 60 K, respectively. The assumption of one-dimensional heat transfer 
is therefore reasonable at the higher strain rate, but horizontal diffusion of heat increases 
temperature significantly at the lower strain rate. If horizontal heat transport were included in the 
thermal models at the lower strain rate, the reaction rate of the phase change would be faster, the 
maximum elevation would be less, and the corresponding upper bound on T c would be lower. 

These considerations favor the formation of Maxwell Montes at horizontal strain rates on the 
order of 10* 15 s* 1 . The implied age of the mountains is of order 50 My or less. Such an age is 
young compared with the average crater retention age of the surface of Venus of about 500 My 
[Phillips et al., 1991]. An important question is whether such a young age is consistent with the 
presence of the 100-km diameter Cleopatra crater [. Basilevsky et al ., 1987] (Figure 1) in Maxwell 
Montes. The age of a single crater, of course, even a large one, is difficult to specify. The rate of 
formation of impact craters larger than 20 km in diameter on Venus has been estimated as 3.3 ±1.8 
x 10*15 km* 2 yr 1 from the statistics of Earth-crossing and Venus-crossing asteroids and an 
assumed set of scaling laws [Shoemaker et al., 1991]. If Maxwell Montes, 1000 km long and 500 
km wide, formed in 50 My, the expected number of craters larger than 20 km in diameter is -0.08. 
The probability that one and only one crater larger than 20 km diameter formed in Maxwell Montes 
in 50 My is thus -8 %. This value indicates that the presence of a large impact crater is not like ly if 
the mountains formed in only 50 My, but the figure is not so low as to reject the hypothesis at high 
confidence. 

Although Danu Montes display compressional deformational features as extensive as the other 
mountain belts of Ishtar Terra, their maximum relief with respect to Lakshmi Planum is as little as 
1 km (Figure 1). The presence of magmatic features within Danu Montes [Solomon etal., 1991, 
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1992; Head et al., 1991; Kaula ei al., 1992] indicates that temperature beneath the mountains 
exceeded the solidus of crustal or shallow mantle material at least locally. Given such evidence for 
relatively high temperature in the lower lithosphere, the reaction of sufficiendy deep lower crustal 
material has probably gone nearly to completion, which may account for the comparatively modest 
topographic relief of these mountains. It is also noteworthy that even a small amount of melting 
can greatly enhance grain-boundary diffusion [Condit et al ., 1985; Watson, 1991]. Assessing the 
cause of higher temperatures beneath Danu Montes requires more detailed thermal models than the 
simple one-dimensional ones considered here. 

CONCLUSIONS 

Taking into account the temperature-dependent reaction rate of the gabbro - eclogite phase 
transition, horizontal strain rates of 10 45 and 10 46 s 4 result in significant differences in the 
maximum elevation of mountains on Venus, not only because of the difference in the formation 
time for relief, but also because of the difference in the thermal regime from advection-dominated 
to crustal-heat-production dominated. For a strain rate of 10 45 s 4 , the observed maximum 
elevation of Maxwell Montes can be explained as a consequence of local isostasy and 
disequilibrium phase boundary depth for a wide range of physical parameters. For the lesser 
horizontal strain rate of 10 46 s* 1 , only limited parameter values for thermal models are allowed, 
and this limitation becomes more difficult to satisfy when a smaller density contrast between the 
crust and mantle is assumed or when horizontal heat transport is considered. We favor the 
formation of Maxwell Montes at a strain rate of lO" 15 s* 1 or greater, although the likelihood that 
Cleopatra crater is younger than the implied age for Maxwell Montes (~50 My) is low. Magmatism 
within Danu Montes indicates higher temperatures in the crust or shallow mantle beneath those 
mountains than beneath the other mountains. The modest elevation of Danu Montes relative to the 
other mountain belts of Ishtar Terra can be understood as a consequence of near completion of the 


phase transition. 
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FIGURE CAPTIONS 

Figure 1. Topographic map of western Ishtar Terra. Contour map of Magellan altimetry data 
gridded at 0.5° intervals of latitude and longitude; 1-km contour interval; polar stereographic 
projection. Cross symbol shows the location of Cleopatra crater (65.9°N 7 .0°E). 

Figure 2. Diffusion constants for ionic species in olivine, pyroxene, and garnet Solid lines show 
diffusion constants that were experimentally determined (1 is for Fe in garnet [Chakraborty and 
Ganguly , 1991]; 2 is for Mg in garnet [Chakraborty and Ganguly, 1991]; 3 is for Fe in olivine 
[Morioka and Nagasawa, 1991]; 4 is for Mg in olivine [Morioka and Nagasawa, 1991]; 5 is for 
Mg in clinopyroxene [Cygan and Lasaga, 1985]; and 6 is for Ca in olivine [Morioka and 
Nagasawa, 1991]). The shaded area shows the range of estimates for Dai,Opx. including 
maximum and minimum (7) values [Smith and Barron, 1991]. Cross symbols show the values 
determined at a single temperature (8 is for Ca in garnet [Chakraborty and Ganguly, 1991]; 9 is an 
upper bound on the diffusion constant of A1 in clinopyroxene [Freer et al., 1988]; 10 and 1 1 are 
values for Al in clinopyroxene from Sautter et al. [1988] and Sautter and Harte [1990], 
respectively). The constants for grain-boundary diffusion of O in forsterite [Condit et al., 1985] 
are converted to equivalent volume diffusion by multiplying by the ratio of grain boundary width to 
grain diameter (12); the grain boundary width was taken to be 3 pm [Condit et al., 1985], and the 
grain diameters were taken to be 1-10 mm. 

Figure 3. Thermal evolution of crust (solid) and mantle (dashed) thickened at a uniform horizontal 
strain rate of (a) 10 45 and (b) 10- 16 s 4 . The phase boundary of the gabbro (G), garnet granulite 
(GG), and eclogite (E) assemblages and the solidus (S) [I to and Kennedy, 1971] are shown by 
dotted lines. The models shown correspond to models (a) 1 and (b) 4 and 5 in Table 1. 

Figure 4. Temporal variation of elevation for the thermal models in Table 1. 


20 


TABLE 1. Parameters for Thermal Models 



Model 

r, mm 

i, s- 1 

D 

T e ,K 

T W ,K 

— 

1 

10 

10-15 

DAl,Opx 

1050 

1321 


2 

1 

10-15 

^Al,Opx 

1120 

1461 


3 

1 

10-15 

DMg,Gt 

1130 

1481 

- 

4 

10 

10-16 

l^ALOpx 

1120 

1461 


5 

1 

10-16 

DAl,Opx 

1120 

1461 

1-^ 

6 

10 

10-16 

Dmr.Gi 

1030 

1281 


m 


m 


a 





Diffusion Constant, m 2 /s 


Fig. 2 
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